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