]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowEvent.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing eff c++ warnings
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowEvent.cxx
CommitLineData
30a892e3 1//////////////////////////////////////////////////////////////////////
2//
3// $Id$
4//
5// Author: Emanuele Simili
6//
7//////////////////////////////////////////////////////////////////////
8//
9//_____________________________________________________________
10//
11// Description:
12// AliFlowEvent is the basic class to perform flow study in ALICE.
13// The AliFlowEvent object contains data mebers wich summarize the ESD
14// informations that are most useful for flow study, together with a
15// collection of tracks (AliFlowTrackCollection) and reconstructed v0s.
16// The class also implements member functions to calculate Multiplicity,
17// Mean P and Pt, Q (the event plane vector), Psi (the event plane angle),
18// normalized q (Q/sqrt(Mult) and Q/sqrt(sum of weights^2)), moreover,
19// functions to separate random or eta sub-events, to fill the bayesian vector
20// of particle abundance (for bayesian p.id calculation), to plug-in weights
21// (to make the tracks' distribution isotropic and to enhance the resolution),
22// and to select the tracks that are used in the event plane calculation (for
23// each harmonic, selection, and subevent.
24// AliFlowEvent supports the standard analysis as well as the Cumulant
25// Analysis: method to extract the generating functions for the old and new
26// cumulant methods is also there.
27// The Flow Event structure (AliFlowEvent, AliFlowTrack, AliFlowV0 and
28// AliFlowSelection) is independent from AliRoot, i.e. once the FlowEvents
29// are filled from the AliESD or the MC KineTree (see AliFlowInterface: here
30// AliRoot is needed), the flow analysis itself can be performed just under
31// ROOT (see AliFlowAnalysisMaker).
32//
33// AliFlowEvent is adapted from the original class StFlowEvent, succesfully
34// employed to study flow in the STAR experiment at RICH.
35// Original Authors: Raimond Snellings & Art Poskanzer
36//
37
38#include "AliFlowEvent.h"
92016a03 39#include "AliFlowTrack.h"
40#include "AliFlowV0.h"
41#include "AliFlowSelection.h"
30a892e3 42#include "AliFlowConstants.h"
92016a03 43
30a892e3 44#include "TVector.h"
45#include "TVector2.h"
46#include "TVector3.h"
92016a03 47#include "TClonesArray.h"
48#include "TMath.h"
49#include "TRandom.h"
50#include "TString.h"
30a892e3 51#include <iostream>
52using namespace std; //required for resolving the 'cout' symbol
53
54// - Flags & Sets
92016a03 55Bool_t AliFlowEvent::fgPtWgt = kFALSE ; // gives pT-based weights
56Bool_t AliFlowEvent::fgEtaWgt = kFALSE ; // gives eta-based weights
4e566f2f 57Bool_t AliFlowEvent::fgOnePhiWgt = kTRUE ; // kTRUE --> ENABLEs SINGLE WEIGHT ISTOGRAM , kFALSE --> ENABLEs 3 WEIGHT ISTOGRAMS
92016a03 58Bool_t AliFlowEvent::fgNoWgt = kFALSE ; // No Weight is used
30a892e3 59// - Eta Sub-Events (later used to calculate the resolution)
4e566f2f 60Bool_t AliFlowEvent::fgEtaSubs = kFALSE ; // makes eta subevents
61Bool_t AliFlowEvent::fgCustomRespFunc = kFALSE ; // custom "detector response function" is used for P.Id (see AliFlowTrack)
30a892e3 62
f847c9e2 63ClassImp(AliFlowEvent)
30a892e3 64//-----------------------------------------------------------
4e566f2f 65AliFlowEvent::AliFlowEvent(Int_t lenght):
d8b98e74 66 fEventID(0), fRunID(0), fOrigMult(0), fL0Trigger(0), fZDCpart(0), fCentrality(-1), fDone(kFALSE),
67 fTrackCollection(0x0), fV0Collection(0x0)
30a892e3 68{
69 // Default constructor: initializes the ObjArray of FlowTracks and FlowV0s,
70 // cleans the internal variables, sets all the weights to 1, sets default flags.
71
30a892e3 72 for(int zz=0;zz<3;zz++) { fZDCenergy[zz] = 0. ; }
92016a03 73 for(int i=0;i<AliFlowConstants::kHars;i++) { fExtPsi[i] = 0. ; fExtRes[i] = 0.; }
30a892e3 74
75 // Make a new track collection
92016a03 76 fTrackCollection = new TClonesArray("AliFlowTrack",lenght) ; fTrackCollection->BypassStreamer(kTRUE) ;
77 fV0Collection = new TClonesArray("AliFlowV0",lenght) ; fV0Collection->BypassStreamer(kTRUE) ;
30a892e3 78
79 // Set Weights Arrays to 1 (default)
327288af 80 for(int nS=0;nS<AliFlowConstants::kSels;nS++)
30a892e3 81 {
327288af 82 for(int nH=0;nH<AliFlowConstants::kHars;nH++)
30a892e3 83 {
327288af 84 for(int nP=0;nP<AliFlowConstants::kPhiBins;nP++)
30a892e3 85 {
86 // enable this with: SetOnePhiWgt()
87 fPhiWgt[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgt[nS][nH][nP] << endl ;
88 // enable these with: SetFirstLastPhiWgt()
89 fPhiWgtPlus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtPlus[nS][nH][nP] << endl ;
90 fPhiWgtMinus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtMinus[nS][nH][nP] << endl ;
91 fPhiWgtCross[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtCross[nS][nH][nP] << endl ;
92 }
93 }
94 }
327288af 95 //for(int nH=0;nH<AliFlowConstants::kHars;nH++) { fExtPsi[nH] = 0. ; fExtRes[nH] = 0. ; }
30a892e3 96
327288af 97 // The Expected particles abundance is taken directly from AliFlowConstants::fgBayesian[] (see Bayesian P.Id.)
30a892e3 98
30a892e3 99}
100//-----------------------------------------------------------
101AliFlowEvent::~AliFlowEvent()
102{
103 // Default distructor: deletes the ObjArrays.
104
105 fTrackCollection->Delete() ; delete fTrackCollection ;
106 fV0Collection->Delete() ; delete fV0Collection ;
107}
108//-------------------------------------------------------------
109
110//-------------------------------------------------------------
111Double_t AliFlowEvent::PhiWeightRaw(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
112{
113 // Weight for making the event plane isotropic in the lab.
114
115 float phi = pFlowTrack->Phi() ; if(phi < 0.) { phi += 2*TMath::Pi() ; }
116 Double_t eta = (Double_t)pFlowTrack->Eta() ;
327288af 117 int n = (int)((phi/(2*TMath::Pi()))*AliFlowConstants::kPhiBins);
30a892e3 118
119 Double_t phiWgt = 1. ;
120 if(OnePhiWgt())
121 {
122 phiWgt = (Double_t)fPhiWgt[selN][harN][n]; //cout << "Selection " << selN << " ; Harmonic " << harN << " ; PhiBin " << n << " - Wgt = " << phiWgt << endl ;
123 }
124 else if(FirstLastPhiWgt())
125 {
126 float zFirstPoint = pFlowTrack->ZFirstPoint(); // float zLastPoint = pFlowTrack->ZLastPoint();
127
128 if (zFirstPoint > 0. && eta > 0.) { phiWgt = (Double_t)fPhiWgtPlus[selN][harN][n] ; }
129 else if(zFirstPoint < 0. && eta < 0.) { phiWgt = (Double_t)fPhiWgtMinus[selN][harN][n] ; }
130 else { phiWgt = (Double_t)fPhiWgtCross[selN][harN][n] ; }
131 }
132
30a892e3 133 return phiWgt ;
134}
135//-------------------------------------------------------------
136Double_t AliFlowEvent::Weight(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
137{
138 // Weight for enhancing the resolution (eta gives sign +/- for Odd Harmonics)
139
327288af 140 if(selN>AliFlowConstants::kSels) { selN = 0 ; }
30a892e3 141 bool oddHar = (harN+1) % 2 ;
142 Double_t phiWgt = 1. ;
143 if(PtWgt())
144 {
145 Double_t pt = (Double_t)pFlowTrack->Pt();
327288af 146 if(pt < AliFlowConstants::fgPtWgtSaturation) { phiWgt *= pt ; }
147 else { phiWgt *= AliFlowConstants::fgPtWgtSaturation ; } // pt weighting going constant
30a892e3 148 }
149 Double_t eta = (Double_t)pFlowTrack->Eta();
150 Double_t etaAbs = TMath::Abs(eta);
151 if(EtaWgt() && oddHar) { phiWgt *= etaAbs ; }
152 if(oddHar && eta < 0.) { phiWgt *= -1. ; }
153
154 return phiWgt ;
155}
156//-------------------------------------------------------------
157Double_t AliFlowEvent::PhiWeight(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
158{
159 // Weight for making the event plane isotropic in the lab and enhancing the resolution
92016a03 160 // (it simply rerurns PhiWeightRaw() * Weight()). If fgNoWgt = kTRUE, returns +/-1 ,
30a892e3 161 // basing on Sign(eta), for odd harmonics .
162
92016a03 163 if(fgNoWgt) // no weights (but +/- according to eta)
30a892e3 164 {
165 bool oddHar = (harN+1) % 2 ;
166 if(oddHar) { return TMath::Sign((Double_t)1.,(Double_t)pFlowTrack->Eta()) ; }
167 else { return (Double_t)1. ; }
168 }
169 Double_t phiWgtRaw = PhiWeightRaw(selN, harN, pFlowTrack);
170 Double_t weight = Weight(selN, harN, pFlowTrack);
327288af 171 if(AliFlowConstants::fgDebug) { cout << "[PhiWeight]: phiWgtRaw = " << phiWgtRaw << " , weight = " << weight << " , eta = " << pFlowTrack->Eta() << endl ; }
30a892e3 172
173 return phiWgtRaw * weight;
174}
175//-------------------------------------------------------------
92016a03 176Int_t AliFlowEvent::Mult(AliFlowSelection* pFlowSelect) const
30a892e3 177{
178 // Multiplicity of tracks in the specified Selection
179
180 if(fDone)
181 {
182 int sub = pFlowSelect->Sub() ;
183 if(sub<0) { return fMult[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
184 else { return fMultSub[sub][pFlowSelect->Sel()][pFlowSelect->Har()] ; }
185 }
186 // -
187 Int_t mult = 0;
188 Int_t itr ;
189 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
190 {
191 AliFlowTrack* pFlowTrack ;
192 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
193 if(pFlowSelect->Select(pFlowTrack)) { mult++ ; }
194 }
195 return mult;
196}
197//-------------------------------------------------------------
92016a03 198Float_t AliFlowEvent::MeanPt(AliFlowSelection* pFlowSelect) const
30a892e3 199{
200 // Mean pt of tracks in the specified Selection
201
202 Double_t meanPt = 0. ;
203 Float_t sumPt = 0. ;
204 UInt_t mult = 0 ;
205 Int_t itr ;
206 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
207 {
208 AliFlowTrack* pFlowTrack ;
209 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
210 if (pFlowSelect->Select(pFlowTrack))
211 {
212 sumPt += pFlowTrack->Pt();
213 mult++;
214 }
215 }
216 if(mult) { meanPt = sumPt/(float)mult ; }
217
218 return meanPt ;
219}
220//-------------------------------------------------------------
92016a03 221TVector2 AliFlowEvent::Q(AliFlowSelection* pFlowSelect) const
30a892e3 222{
223 // Event plane vector for the specified Selection
224
225 if(fDone)
226 {
227 int sub = pFlowSelect->Sub() ;
228 if(sub<0) { return fQ[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
229 else { return fQSub[sub][pFlowSelect->Sel()][pFlowSelect->Har()] ; }
230 }
231 // -
232 TVector2 mQ ;
233 Double_t mQx=0. , mQy=0. ;
234 int selN = pFlowSelect->Sel() ;
235 int harN = pFlowSelect->Har() ;
236 double order = (double)(harN + 1) ;
237
238 Int_t itr ;
239 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
240 {
241 AliFlowTrack* pFlowTrack ;
242 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
243 if(pFlowSelect->Select(pFlowTrack))
244 {
245 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
246 float phi = pFlowTrack->Phi();
247 mQx += phiWgt * cos(phi * order) ;
248 mQy += phiWgt * sin(phi * order) ;
327288af 249 if(AliFlowConstants::fgDebug) { cout << itr << " phi = " << phi << " , wgt = " << phiWgt << endl ; }
30a892e3 250 }
251 }
252 mQ.Set(mQx, mQy);
253
254 return mQ;
255}
256//-------------------------------------------------------------
92016a03 257Float_t AliFlowEvent::Psi(AliFlowSelection* pFlowSelect) const
30a892e3 258{
259 // Event plane angle for the specified Selection
260
261 int harN = pFlowSelect->Har() ;
262 float order = (float)(harN + 1) ;
263 Float_t psi = 0. ;
264
265 TVector2 mQ = Q(pFlowSelect);
266 if(mQ.Mod()) // if vector is not 0
267 {
268 psi= mQ.Phi() / order ;
269 if (psi < 0.) { psi += 2*TMath::Pi() / order ; }
270 }
271 return psi;
272}
273//-------------------------------------------------------------
92016a03 274TVector2 AliFlowEvent::NormQ(AliFlowSelection* pFlowSelect) const
30a892e3 275{
276 // Normalized Q = Q/sqrt(sum of weights^2) for the specified Selection
277
278 if(fDone)
279 {
280 TVector2 mQ = fQ[pFlowSelect->Sel()][pFlowSelect->Har()] ;
92016a03 281 double sumOfWeightSqr = fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ;
282 if(sumOfWeightSqr) { mQ /= TMath::Sqrt(sumOfWeightSqr) ; }
30a892e3 283 else { mQ.Set(0.,0.) ; }
284 return mQ ;
285 }
286 // -
287 TVector2 mQ ;
288 Double_t mQx=0. , mQy=0. ;
92016a03 289 double sumOfWeightSqr = 0 ;
30a892e3 290 int selN = pFlowSelect->Sel() ;
291 int harN = pFlowSelect->Har() ;
292 double order = (double)(harN + 1) ;
293 Int_t itr ;
294 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
295 {
296 AliFlowTrack* pFlowTrack ;
297 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
298 if (pFlowSelect->Select(pFlowTrack))
299 {
300 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
92016a03 301 sumOfWeightSqr += phiWgt*phiWgt;
30a892e3 302
303 float phi = pFlowTrack->Phi();
304 mQx += phiWgt * cos(phi * order);
305 mQy += phiWgt * sin(phi * order);
306 }
307 }
92016a03 308 if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
30a892e3 309 else { mQ.Set(0.,0.); }
310
311 return mQ;
312}
313//-------------------------------------------------------------
92016a03 314Float_t AliFlowEvent::OldQ(AliFlowSelection* pFlowSelect) const
30a892e3 315{
316 // Magnitude of normalized Q vector (without pt or eta weighting) for the specified Selection
317
318 TVector2 mQ ;
319 Double_t mQx = 0. , mQy = 0. ;
320 int selN = pFlowSelect->Sel() ;
321 int harN = pFlowSelect->Har() ;
322 double order = (double)(harN + 1) ;
92016a03 323 double sumOfWeightSqr = 0 ;
30a892e3 324
325 Int_t itr ;
326 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
327 {
328 AliFlowTrack* pFlowTrack ;
329 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
330 if(pFlowSelect->Select(pFlowTrack))
331 {
332 double phiWgt = PhiWeightRaw(selN, harN, pFlowTrack); // Raw
92016a03 333 sumOfWeightSqr += phiWgt*phiWgt;
30a892e3 334 float phi = pFlowTrack->Phi();
335 mQx += phiWgt * cos(phi * order);
336 mQy += phiWgt * sin(phi * order);
337 }
338 }
92016a03 339 if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
30a892e3 340 else { mQ.Set(0.,0.); }
341
342 return mQ.Mod();
343}
344//-----------------------------------------------------------------------
92016a03 345Double_t AliFlowEvent::NewG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
30a892e3 346{
347 // Generating function for the new cumulant method. Eq. 3 in the Practical Guide
348
349 int selN = pFlowSelect->Sel();
350 int harN = pFlowSelect->Har();
351 double order = (double)(harN + 1);
352
353 double mult = (double)Mult(pFlowSelect);
354 Double_t theG = 1.;
355
356 Int_t itr ;
357 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
358 {
359 AliFlowTrack* pFlowTrack ;
360 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
361 if (pFlowSelect->Select(pFlowTrack))
362 {
363 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
364 float phi = pFlowTrack->Phi();
365 theG *= (1. + (phiWgt/mult) * (2.* Zx * cos(phi * order) + 2.* Zy * sin(phi * order) ) );
366 }
367 }
368 return theG;
369}
370//-----------------------------------------------------------------------
92016a03 371Double_t AliFlowEvent::OldG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
30a892e3 372{
373 // Generating function for the old cumulant method (if expanded in Taylor
92016a03 374 // series, one recovers NewG() in new new cumulant method)
30a892e3 375
92016a03 376 TVector2 normQ = NormQ(pFlowSelect) ;
30a892e3 377
378 return exp(2*Zx*normQ.X() + 2*Zy*normQ.Y());
379}
380//-----------------------------------------------------------------------
92016a03 381Double_t AliFlowEvent::SumWeightSquare(AliFlowSelection* pFlowSelect) const
30a892e3 382{
383 // Return sum of weights^2 for the specified Selection (used for normalization)
384
92016a03 385 if(fDone) { return fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
386
30a892e3 387 int selN = pFlowSelect->Sel();
388 int harN = pFlowSelect->Har();
92016a03 389 Double_t sumOfWeightSqr = 0;
30a892e3 390 Int_t itr ;
391 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
392 {
393 AliFlowTrack* pFlowTrack ;
394 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
395 if (pFlowSelect->Select(pFlowTrack))
396 {
397 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
92016a03 398 sumOfWeightSqr += phiWgt*phiWgt;
30a892e3 399 }
400 }
92016a03 401 if(sumOfWeightSqr < 0.) { return Mult(pFlowSelect) ; }
30a892e3 402
92016a03 403 return sumOfWeightSqr;
30a892e3 404}
405//-------------------------------------------------------------
92016a03 406Double_t AliFlowEvent::WgtMultQ4(AliFlowSelection* pFlowSelect) const
30a892e3 407{
408 // Used only for the old cumulant method, for getting q4 when weight is on.
409 // Replace multiplicity in Eq.(74b) by this quantity when weight is on.
410 // This is derived based on (A4) in the old cumulant paper.
411
412 int selN = pFlowSelect->Sel();
413 int harN = pFlowSelect->Har();
414 double theMult = 0.;
415 double theMeanWj4 = 0.;
416 double theMeanWj2 = 0.;
417 double theSumOfWgtSqr = 0;
418 double phiWgtSq;
419
420 Int_t itr ;
421 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
422 {
423 AliFlowTrack* pFlowTrack ;
424 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
425 if (pFlowSelect->Select(pFlowTrack))
426 {
427 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
428 phiWgtSq = phiWgt*phiWgt;
429 theSumOfWgtSqr += phiWgtSq;
430 theMeanWj4 += phiWgtSq*phiWgtSq;
431 theMult += 1.;
432 }
433 }
434 if (theMult <= 0.) return theMult;
435
436 theMeanWj4 /= theMult;
437 theMeanWj2 = theSumOfWgtSqr / theMult;
438
439 return (theSumOfWgtSqr*theSumOfWgtSqr)/(theMult*(-theMeanWj4+2*theMeanWj2*theMeanWj2));
440}
441//-------------------------------------------------------------
92016a03 442Double_t AliFlowEvent::WgtMultQ6(AliFlowSelection* pFlowSelect) const
30a892e3 443{
444 // Used only for the old cumulant method. For getting q6 when weight is on.
445 // Replace multiplicity in Eq.(74c) by this quantity when weight is on.
446 // This is derived based on (A4) in the old cumulant paper.
447
448 int selN = pFlowSelect->Sel();
449 int harN = pFlowSelect->Har();
450 double theMult = 0.;
451 double theMeanWj6 = 0.;
452 double theMeanWj4 = 0.;
453 double theMeanWj2 = 0.;
454 double theSumOfWgtSqr = 0;
455 double phiWgtSq;
456
457 Int_t itr ;
458 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
459 {
460 AliFlowTrack* pFlowTrack ;
461 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
462 if (pFlowSelect->Select(pFlowTrack))
463 {
464 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
465 phiWgtSq = phiWgt*phiWgt;
466 theSumOfWgtSqr += phiWgtSq;
467 theMeanWj4 += phiWgtSq*phiWgtSq;
468 theMeanWj6 += phiWgtSq*phiWgtSq*phiWgtSq;
469 theMult += 1.;
470 }
471 }
472 if (theMult <= 0.) return theMult*theMult;
473
474 theMeanWj6 /= theMult;
475 theMeanWj4 /= theMult;
476 theMeanWj2 = theSumOfWgtSqr / theMult;
477
478 return 4.*(theSumOfWgtSqr*theSumOfWgtSqr*theSumOfWgtSqr)/(theMult*(theMeanWj6-9.*theMeanWj2*theMeanWj4+12.*theMeanWj2*theMeanWj2*theMeanWj2));
479}
480//-------------------------------------------------------------
481void AliFlowEvent::SetSelections(AliFlowSelection* pFlowSelect)
482{
483 // Sets the selection of tracks used in the Reaction Plane calculation
484 // for the specific Harmonic and Selection - this does not cut trow away
485 // tracks from the event, neither exclude them from flow study. See also
486 // the AliFlowSelection class.
487 // Strategy of Selection:
488 // For the specific Harmonic and Selection, IF cuts
489 // are defined (such that low<high) and IF the track satisfies them, THEN
490 // the respective track's flag (bool AliFlowTrack::mSelection[har][sel])
491 // is set kTRUE so that the track, from now on, will be included in the
492 // R.P. determination for that selection.
493 // If NO cuts are defined -> ALL Flags are setted kTRUE (all particle
494 // used for all the Reaction Planes -> no difference in Psi[har][sel]).
495 // -------------------------------------------------------------------
496 // The first selection (all harmonics) is set kTRUE : no conditions.
497
498 Int_t itr ;
499 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
500 {
501 AliFlowTrack* pFlowTrack ;
502 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
503 pFlowTrack->ResetSelection() ; // this re-sets all the mSelection flags to 0
504
505 // * this sets all the selection n.[0] flag kTRUE (all harmonics) *
327288af 506 for(int harN=0;harN<AliFlowConstants::kHars;harN++) { pFlowTrack->SetSelect(harN,0) ; }
30a892e3 507
508 // Track need to be Constrainable
509 if(pFlowSelect->ConstrainCut() && !pFlowTrack->IsConstrainable()) continue ;
510
511 // PID - gets rid of the track with AliFlowTrack::Pid() != AliFlowSelection::Pid() (if there)
512 if(pFlowSelect->Pid()[0] != '\0')
513 {
514 if(strstr(pFlowSelect->Pid(), "h")!=0)
515 {
516 int charge = pFlowTrack->Charge() ;
517 if(strcmp("h+",pFlowSelect->Pid())==0 && charge != 1) continue;
518 if(strcmp("h-",pFlowSelect->Pid())==0 && charge != -1) continue;
519 }
520 else
521 {
522 Char_t pid[10];
523 strcpy(pid, pFlowTrack->Pid());
524 if(strstr(pid, pFlowSelect->Pid())==0) continue;
525 }
526 }
527 double eta = (double)(pFlowTrack->Eta());
92016a03 528 float pt = pFlowTrack->Pt();
30a892e3 529 float gDca = pFlowTrack->Dca() ;
530
531 // Global DCA - gets rid of the track with DCA outside the range
532 if((pFlowSelect->DcaGlobalCutHi()>pFlowSelect->DcaGlobalCutLo()) && (gDca<pFlowSelect->DcaGlobalCutLo() || gDca>pFlowSelect->DcaGlobalCutHi())) continue ;
533
534 // Pt & Eta - this is done differently for different Harmonic & Selection
327288af 535 for(int selN = 1; selN < AliFlowConstants::kSels; selN++) // not even consider the 0th selection (no cut applied there)
30a892e3 536 {
537 // min. TPC hits required
538 if(pFlowSelect->NhitsCut(selN) && (pFlowTrack->FitPtsTPC()<pFlowSelect->NhitsCut(selN))) continue ;
539
327288af 540 for(int harN = 0; harN < AliFlowConstants::kHars; harN++)
30a892e3 541 {
542 // Eta - gets rid of the track with Eta outside the range
543 if((pFlowSelect->EtaCutHi(harN%2,selN)>pFlowSelect->EtaCutLo(harN%2,selN)) && (TMath::Abs(eta)<pFlowSelect->EtaCutLo(harN%2,selN) || TMath::Abs(eta)>pFlowSelect->EtaCutHi(harN%2,selN))) continue ;
544 // Pt - gets rid of the track with Pt outside the range
92016a03 545 if((pFlowSelect->PtCutHi(harN%2,selN)>pFlowSelect->PtCutLo(harN%2,selN)) && (pt<pFlowSelect->PtCutLo(harN%2,selN) || pt>pFlowSelect->PtCutHi(harN%2,selN))) continue ;
30a892e3 546
547 pFlowTrack->SetSelect(harN, selN) ; // if cuts defined (low<high) && track is in the range -> Set [har][sel] Flag ON
548
327288af 549 if(AliFlowConstants::fgDebug)
30a892e3 550 {
92016a03 551 if(harN==1)
552 {
553 cout << " n. " << itr << " , pFlowTrack->PrintSelection() = " ; pFlowTrack->PrintSelection() ;
554 if(pFlowSelect->Pid()[0] != '\0') { cout << " track: pid " << pFlowTrack->Pid() << " = "<< pFlowSelect->Pid() << endl ; }
555 cout << " track: dca " << pFlowSelect->DcaGlobalCutLo() << " < " << gDca << " < " << pFlowSelect->DcaGlobalCutHi() << endl ;
556 cout << " track: eta " << pFlowSelect->EtaCutLo(harN,selN) << " < |" << eta << "| < " << pFlowSelect->EtaCutHi(harN,selN) << endl ;
557 cout << " track: pT " << pFlowSelect->PtCutLo(harN,selN) << " < " << pt << " < " << pFlowSelect->PtCutHi(harN,selN) << endl ;
558 cout << " pFlowTrack->PrintSelection() = " ; pFlowTrack->PrintSelection() ;
559 }
560 // cout << " selN " << selN << " , harN " << harN%2 << " - si" << endl ;
30a892e3 561 }
562 }
563 }
564 }
565}
566//-------------------------------------------------------------
567void AliFlowEvent::SetPids()
568{
4e566f2f 569 // Re-sets the tracks P.id. (using the current AliFlowConstants::fgBayesian[] array).
570 // To use the Raw P.Id (just the detector response function), set fgBayesian[] = {1,1,1,1,1,1}.
30a892e3 571
92016a03 572 const Int_t kCode[] = {11,13,211,321,2212,10010020} ;
30a892e3 573 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
574 {
575 AliFlowTrack* pFlowTrack ;
576 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
4e566f2f 577
578 Float_t bayPid[AliFlowConstants::kPid] ;
579 if(fgCustomRespFunc) { pFlowTrack->PidProbsC(bayPid) ; }
580 else { pFlowTrack->PidProbs(bayPid) ; }
581
582 Int_t maxN = 2 ; // if No id. -> then is a Pi
583 Float_t pidMax = bayPid[2] ; // (if all equal, Pi probability was the first)
327288af 584 for(Int_t nP=0;nP<AliFlowConstants::kPid;nP++)
30a892e3 585 {
92016a03 586 if(bayPid[nP]>pidMax) { maxN = nP ; pidMax = bayPid[nP] ; }
30a892e3 587 }
92016a03 588 Int_t pdgCode = TMath::Sign(kCode[maxN],pFlowTrack->Charge()) ;
4e566f2f 589 pFlowTrack->SetMostLikelihoodPID(pdgCode) ;
30a892e3 590 }
591}
592//-------------------------------------------------------------
92016a03 593void AliFlowEvent::MakeSubEvents() const
30a892e3 594{
595 // Make random or eta sub-events
596
92016a03 597 if(fgEtaSubs) { MakeEtaSubEvents() ; }
30a892e3 598 else { MakeRndSubEvents() ; }
599}
600//-------------------------------------------------------------
92016a03 601void AliFlowEvent::MakeRndSubEvents() const
30a892e3 602{
603 // Make random subevents
604
327288af 605 int eventMult[AliFlowConstants::kHars][AliFlowConstants::kSels] = {{0}};
30a892e3 606 int harN, selN, subN = 0;
607
608 // loop to count the total number of tracks for each selection
609 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
610 {
611 AliFlowTrack* pFlowTrack ;
612 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
327288af 613 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
30a892e3 614 {
327288af 615 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
30a892e3 616 {
617 if(pFlowTrack->Select(harN, selN)) { eventMult[harN][selN]++ ; }
618 }
619 }
620 }
621 // loop to set the SubEvent member
327288af 622 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
30a892e3 623 {
327288af 624 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
30a892e3 625 {
327288af 626 int subEventMult = eventMult[harN][selN] / AliFlowConstants::kSubs;
30a892e3 627 if (subEventMult)
628 {
629 subN = 0;
630 int countN = 0;
631 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
632 {
633 AliFlowTrack* pFlowTrack ;
634 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
635 if(pFlowTrack->Select(harN, selN))
636 {
637 pFlowTrack->SetSubevent(harN, selN, subN);
638 countN++;
639 if (countN % subEventMult == 0.) { subN++ ; }
640 }
641 }
642 }
643 }
644 }
645 return ;
646}
647//-------------------------------------------------------------
92016a03 648void AliFlowEvent::MakeEtaSubEvents() const
30a892e3 649{
650 // Make subevents for positive and negative eta
92016a03 651 // (when done, fgEtaSubs flag setted to kTRUE).
30a892e3 652
653 int harN, selN = 0;
654 // loop to set the SubEvent member
327288af 655 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
30a892e3 656 {
327288af 657 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
30a892e3 658 {
659 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
660 {
661 AliFlowTrack* pFlowTrack ;
662 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
663 if(pFlowTrack->Select(harN, selN))
664 {
665 float eta = pFlowTrack->Eta();
666 if (eta > 0.) { pFlowTrack->SetSubevent(harN, selN, 0) ; }
667 else { pFlowTrack->SetSubevent(harN, selN, 1) ; }
668 }
669 }
670 }
671 }
672}
673//-------------------------------------------------------------
674void AliFlowEvent::RandomShuffle()
675{
676 // Randomly re-shuffles the tracks in the array; if a track is not
677 // primary, the reference carried by the reconstructed mother (in
678 // the v0 array) is changed accordingly.
679
92016a03 680 return ; // ...at the moment is disabled ... //
681
682// Int_t tot = 0 ;
683// UInt_t imax = fTrackCollection->GetEntries() ;
684// TRandom* rnd = new TRandom(0) ; // SetSeed(0) ;
685// TClonesArray* newTrackCollection = new TClonesArray("AliFlowTrack",imax) ;
686//
687// // re-arranges the ObjArray (TrackCollection())
688// for(UInt_t itr=0;itr<imax;itr++)
689// {
690// AliFlowTrack* pFlowTrack ;
691// pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
692//
693// UInt_t rndNumber = rnd->Integer(imax) ;
694// Bool_t put = kFALSE ;
695// while(!put)
696// {
697// if(!newTrackCollection->AddrAt(rndNumber))
698// {
699// ... new(newTrackCollection[rndNumber]) AliFlowTrack(*pFlowTrack) ;
700// put = kTRUE ; tot++ ;
701// if(AliFlowConstants::fgDebug) { cout << " " << itr << " --> " << rndNumber << endl ; }
702// }
703// else
704// {
705// rndNumber++ ; if(rndNumber>=imax) { rndNumber -= imax ; }
706// }
707// }
708// }
709// if(AliFlowConstants::fgDebug) { cout << "* RandomShuffle() : " << tot << "/" << imax << " flow tracks have been shuffled " << endl ; }
710// fTrackCollection = newTrackCollection ;
30a892e3 711}
712//-----------------------------------------------------------------------
713UInt_t AliFlowEvent::Centrality()
714{
715 // Returns the Centrality Class as stored
716
717 if(fCentrality < 0) { SetCentrality() ; }
718 return fCentrality ;
719}
720//-----------------------------------------------------------------------
30a892e3 721void AliFlowEvent::SetCentrality()
722{
723 // Sets the Centrality Classes basing on Multiplicity at mid rapidity,
724 // ... (ZDC information can be added) .
725
726 Float_t* cent ;
727 Int_t tracks = MultEta() ;
728
729 if(RunID() == -1)
730 {
327288af 731 cent = AliFlowConstants::fgCentNorm ;
30a892e3 732 //if centrality classes are not defined, does it now (with CentNorm & MaxMult)
327288af 733 if(cent[AliFlowConstants::kCents-1] <= 1)
30a892e3 734 {
327288af 735 for(Int_t ic=0;ic<AliFlowConstants::kCents;ic++)
30a892e3 736 {
327288af 737 cent[ic] *= AliFlowConstants::fgMaxMult ;
738 if(AliFlowConstants::fgDebug) { cout << "Centrality[" << ic << "] = " << cent[ic] << " . " << endl ; }
30a892e3 739 }
740 }
741 }
742 else if((RunID() != -1) && (CenterOfMassEnergy() == 5500.))
743 {
327288af 744 cent = (Float_t*)AliFlowConstants::fgCent0 ;
30a892e3 745 }
92016a03 746 else // other definition of centrality are possible ...
30a892e3 747 {
327288af 748 cent = (Float_t*)AliFlowConstants::fgCent0 ;
30a892e3 749 }
750 if (tracks < cent[0]) { fCentrality = 0; }
751 else if (tracks < cent[1]) { fCentrality = 1; }
752 else if (tracks < cent[2]) { fCentrality = 2; }
753 else if (tracks < cent[3]) { fCentrality = 3; }
754 else if (tracks < cent[4]) { fCentrality = 4; }
755 else if (tracks < cent[5]) { fCentrality = 5; }
756 else if (tracks < cent[6]) { fCentrality = 6; }
757 else if (tracks < cent[7]) { fCentrality = 7; }
758 else if (tracks < cent[8]) { fCentrality = 8; }
759 else { fCentrality = 9; }
760
327288af 761 if(AliFlowConstants::fgDebug) { cout << " * Centrality Class : " << fCentrality << " . " << endl ; }
30a892e3 762}
763//-----------------------------------------------------------------------
92016a03 764TVector AliFlowEvent::Bayesian() const
30a892e3 765{
92016a03 766 // Returns bayesian array of particle abundances (from AliFlowConstants::)
767
327288af 768 TVector bayes(AliFlowConstants::kPid) ;
769 for(Int_t i=0;i<AliFlowConstants::kPid;i++) { bayes[i] = AliFlowConstants::fgBayesian[i] ; }
30a892e3 770 return bayes ;
771}
772//-----------------------------------------------------------------------
773void AliFlowEvent::PrintFlagList() const
774{
775 // Prints the list of selection cuts ( called in AliFlowInterface::Finish() )
776
777 cout << "#######################################################" << endl;
778 cout << "# Weighting and Striping:" << endl;
779 if(PtWgt())
780 {
781 cout << "# PtWgt = kTRUE " << endl ; // (also for output of PhiWgt file?)
327288af 782 cout << "# PtWgt Saturation = " << AliFlowConstants::fgPtWgtSaturation << endl;
30a892e3 783 }
784 else
785 {
786 cout << "# PtWgt = kFALSE" << endl;
787 }
788 if(EtaWgt())
789 {
790 cout << "# EtaWgt = kTRUE " << endl ; // (also for output of PhiWgt file for odd harmonics?)
791 }
792 else
793 {
794 cout << "# EtaWgt = kFALSE" << endl;
795 }
796 cout << "#######################################################" << endl;
797}
798//-----------------------------------------------------------------------
92016a03 799Int_t AliFlowEvent::MultEta() const
30a892e3 800{
327288af 801 // Returns the multiplicity in the interval |eta|<(AliFlowConstants::fgEetaMid), used
30a892e3 802 // for centrality measurement (see centrality classes in fCentrality) .
803
804 Int_t goodtracks = 0 ;
805 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
806 {
807 AliFlowTrack* pFlowTrack ;
808 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
327288af 809 if((pFlowTrack->Charge()) && (TMath::Abs(pFlowTrack->Eta())<AliFlowConstants::fgEtaMid)) { goodtracks++ ; }
30a892e3 810 }
811 return goodtracks ;
812}
813//-----------------------------------------------------------------------
814Int_t AliFlowEvent::UncorrNegMult(Float_t eta) const
815{
816 // Negative multiplicity in the interval (-eta..eta)
327288af 817 // (default is AliFlowConstants::fgEetaGood = 0.9)
30a892e3 818
819 Int_t negMult = 0 ;
820 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
821 {
822 AliFlowTrack* pFlowTrack ;
823 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
824 if((pFlowTrack->Charge()<0) && (TMath::Abs(pFlowTrack->Eta())<TMath::Abs(eta))) { negMult++ ; }
825 delete pFlowTrack ;
826 }
827 return negMult;
828}
829//-----------------------------------------------------------------------
830Int_t AliFlowEvent::UncorrPosMult(Float_t eta) const
831{
832 // Positive multiplicity in the interval (-eta..eta)
327288af 833 // (default is AliFlowConstants::fgEetaGood = 0.9)
30a892e3 834
51a9b7d8 835 Int_t posMult = 0 ;
30a892e3 836 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
837 {
838 AliFlowTrack* pFlowTrack ;
839 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
840 if((pFlowTrack->Charge()>0) && (TMath::Abs(pFlowTrack->Eta())<TMath::Abs(eta))) { posMult++ ; }
841 delete pFlowTrack ;
842 }
843 return posMult;
844}
845//-----------------------------------------------------------------------
30a892e3 846TVector3 AliFlowEvent::VertexPos() const
847{
92016a03 848 // Returns primary vertex position as a TVector3
849
30a892e3 850 Float_t vertex[3] ;
851 VertexPos(vertex) ;
852 return TVector3(vertex) ;
853}
854//-----------------------------------------------------------------------
30a892e3 855void AliFlowEvent::MakeAll()
856{
92016a03 857 // calculates all quantities in 1 shoot
30a892e3 858
327288af 859 Double_t mQx[AliFlowConstants::kSels][AliFlowConstants::kHars] ;
860 Double_t mQy[AliFlowConstants::kSels][AliFlowConstants::kHars] ;
861 Double_t mQxSub[AliFlowConstants::kSubs][AliFlowConstants::kSels][AliFlowConstants::kHars] ;
862 Double_t mQySub[AliFlowConstants::kSubs][AliFlowConstants::kSels][AliFlowConstants::kHars] ;
30a892e3 863 // -
864 int selN, harN, subN ;
327288af 865 for(selN=0;selN<AliFlowConstants::kSels;selN++)
30a892e3 866 {
327288af 867 for(harN=0;harN<AliFlowConstants::kHars;harN++)
30a892e3 868 {
869 mQx[selN][harN] = 0. ;
870 mQy[selN][harN] = 0. ;
871 fMult[selN][harN] = 0 ;
872 fSumOfWeightSqr[selN][harN] = 0. ;
327288af 873 for(subN=0;subN<AliFlowConstants::kSubs;subN++)
30a892e3 874 {
875 mQxSub[subN][selN][harN] = 0. ;
876 mQySub[subN][selN][harN] = 0. ;
877 fMultSub[subN][selN][harN] = 0 ;
878 }
879 }
880 }
881
882 double order = 0. ;
883 double phiWgt = 0. ;
884 float phi = 0. ;
885 // -
886 int itr ;
887 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
888 {
889 AliFlowTrack* pFlowTrack ;
890 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
891 phi = pFlowTrack->Phi();
327288af 892 for(selN=0;selN<AliFlowConstants::kSels;selN++)
30a892e3 893 {
327288af 894 for(harN=0;harN<AliFlowConstants::kHars;harN++)
30a892e3 895 {
896 order = (double)(harN+1) ;
897 if(pFlowTrack->Select(harN,selN))
898 {
899 phiWgt = PhiWeight(selN,harN,pFlowTrack) ;
900 fSumOfWeightSqr[selN][harN] += phiWgt*phiWgt ;
901 mQx[selN][harN] += phiWgt * cos(phi * order) ;
902 mQy[selN][harN] += phiWgt * sin(phi * order) ;
903 fMult[selN][harN]++ ;
327288af 904 for(subN=0;subN<AliFlowConstants::kSubs;subN++)
30a892e3 905 {
906 if(pFlowTrack->Select(harN,selN,subN))
907 {
908 mQxSub[subN][selN][harN] += phiWgt * cos(phi * order) ;
909 mQySub[subN][selN][harN] += phiWgt * sin(phi * order) ;
910 fMultSub[subN][selN][harN]++ ;
911 }
912 } // sub
913 } // if
914 } // har
915 } // sel
916 } //itr
917
327288af 918 for(selN=0;selN<AliFlowConstants::kSels;selN++)
30a892e3 919 {
327288af 920 for(harN=0;harN<AliFlowConstants::kHars;harN++)
30a892e3 921 {
922 fQ[selN][harN].Set(mQx[selN][harN],mQy[selN][harN]) ;
327288af 923 for(subN=0;subN<AliFlowConstants::kSubs;subN++) { fQSub[subN][selN][harN].Set(mQxSub[subN][selN][harN],mQySub[subN][selN][harN]) ; }
30a892e3 924 }
925 }
926
927 fDone = kTRUE ;
928}
30a892e3 929//-----------------------------------------------------------------------