]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG2/FLOW/AliFlowEvent.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
removing eff c++ warnings - more to come
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowEvent.cxx
... / ...
CommitLineData
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"
39#include "AliFlowTrack.h"
40#include "AliFlowV0.h"
41#include "AliFlowSelection.h"
42#include "AliFlowConstants.h"
43
44#include "TVector.h"
45#include "TVector2.h"
46#include "TVector3.h"
47#include "TClonesArray.h"
48#include "TMath.h"
49#include "TRandom.h"
50#include "TString.h"
51#include <iostream>
52using namespace std; //required for resolving the 'cout' symbol
53
54// - Flags & Sets
55Bool_t AliFlowEvent::fgPtWgt = kFALSE ; // gives pT-based weights
56Bool_t AliFlowEvent::fgEtaWgt = kFALSE ; // gives eta-based weights
57Bool_t AliFlowEvent::fgOnePhiWgt = kTRUE ; // kTRUE --> ENABLEs SINGLE WEIGHT ISTOGRAM , kFALSE --> ENABLEs 3 WEIGHT ISTOGRAMS
58Bool_t AliFlowEvent::fgNoWgt = kFALSE ; // No Weight is used
59// - Eta Sub-Events (later used to calculate the resolution)
60Bool_t AliFlowEvent::fgEtaSubs = kFALSE ; // makes eta subevents
61Bool_t AliFlowEvent::fgCustomRespFunc = kFALSE ; // custom "detector response function" is used for P.Id (see AliFlowTrack)
62
63ClassImp(AliFlowEvent)
64//-----------------------------------------------------------
65AliFlowEvent::AliFlowEvent(Int_t lenght):
66 fEventID(0), fRunID(0), fOrigMult(0), fL0Trigger(0), fZDCpart(0), fCentrality(-1), fDone(kFALSE),
67 fTrackCollection(0x0), fV0Collection(0x0)
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
72 for(int zz=0;zz<3;zz++) { fZDCenergy[zz] = 0. ; }
73 for(int i=0;i<AliFlowConstants::kHars;i++) { fExtPsi[i] = 0. ; fExtRes[i] = 0.; }
74
75 // Make a new track collection
76 fTrackCollection = new TClonesArray("AliFlowTrack",lenght) ; fTrackCollection->BypassStreamer(kTRUE) ;
77 fV0Collection = new TClonesArray("AliFlowV0",lenght) ; fV0Collection->BypassStreamer(kTRUE) ;
78
79 // Set Weights Arrays to 1 (default)
80 for(int nS=0;nS<AliFlowConstants::kSels;nS++)
81 {
82 for(int nH=0;nH<AliFlowConstants::kHars;nH++)
83 {
84 for(int nP=0;nP<AliFlowConstants::kPhiBins;nP++)
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 }
95 //for(int nH=0;nH<AliFlowConstants::kHars;nH++) { fExtPsi[nH] = 0. ; fExtRes[nH] = 0. ; }
96
97 // The Expected particles abundance is taken directly from AliFlowConstants::fgBayesian[] (see Bayesian P.Id.)
98
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() ;
117 int n = (int)((phi/(2*TMath::Pi()))*AliFlowConstants::kPhiBins);
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
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
140 if(selN>AliFlowConstants::kSels) { selN = 0 ; }
141 bool oddHar = (harN+1) % 2 ;
142 Double_t phiWgt = 1. ;
143 if(PtWgt())
144 {
145 Double_t pt = (Double_t)pFlowTrack->Pt();
146 if(pt < AliFlowConstants::fgPtWgtSaturation) { phiWgt *= pt ; }
147 else { phiWgt *= AliFlowConstants::fgPtWgtSaturation ; } // pt weighting going constant
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
160 // (it simply rerurns PhiWeightRaw() * Weight()). If fgNoWgt = kTRUE, returns +/-1 ,
161 // basing on Sign(eta), for odd harmonics .
162
163 if(fgNoWgt) // no weights (but +/- according to eta)
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);
171 if(AliFlowConstants::fgDebug) { cout << "[PhiWeight]: phiWgtRaw = " << phiWgtRaw << " , weight = " << weight << " , eta = " << pFlowTrack->Eta() << endl ; }
172
173 return phiWgtRaw * weight;
174}
175//-------------------------------------------------------------
176Int_t AliFlowEvent::Mult(AliFlowSelection* pFlowSelect) const
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//-------------------------------------------------------------
198Float_t AliFlowEvent::MeanPt(AliFlowSelection* pFlowSelect) const
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//-------------------------------------------------------------
221TVector2 AliFlowEvent::Q(AliFlowSelection* pFlowSelect) const
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) ;
249 if(AliFlowConstants::fgDebug) { cout << itr << " phi = " << phi << " , wgt = " << phiWgt << endl ; }
250 }
251 }
252 mQ.Set(mQx, mQy);
253
254 return mQ;
255}
256//-------------------------------------------------------------
257Float_t AliFlowEvent::Psi(AliFlowSelection* pFlowSelect) const
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//-------------------------------------------------------------
274TVector2 AliFlowEvent::NormQ(AliFlowSelection* pFlowSelect) const
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()] ;
281 double sumOfWeightSqr = fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ;
282 if(sumOfWeightSqr) { mQ /= TMath::Sqrt(sumOfWeightSqr) ; }
283 else { mQ.Set(0.,0.) ; }
284 return mQ ;
285 }
286 // -
287 TVector2 mQ ;
288 Double_t mQx=0. , mQy=0. ;
289 double sumOfWeightSqr = 0 ;
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);
301 sumOfWeightSqr += phiWgt*phiWgt;
302
303 float phi = pFlowTrack->Phi();
304 mQx += phiWgt * cos(phi * order);
305 mQy += phiWgt * sin(phi * order);
306 }
307 }
308 if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
309 else { mQ.Set(0.,0.); }
310
311 return mQ;
312}
313//-------------------------------------------------------------
314Float_t AliFlowEvent::OldQ(AliFlowSelection* pFlowSelect) const
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) ;
323 double sumOfWeightSqr = 0 ;
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
333 sumOfWeightSqr += phiWgt*phiWgt;
334 float phi = pFlowTrack->Phi();
335 mQx += phiWgt * cos(phi * order);
336 mQy += phiWgt * sin(phi * order);
337 }
338 }
339 if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
340 else { mQ.Set(0.,0.); }
341
342 return mQ.Mod();
343}
344//-----------------------------------------------------------------------
345Double_t AliFlowEvent::NewG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
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//-----------------------------------------------------------------------
371Double_t AliFlowEvent::OldG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
372{
373 // Generating function for the old cumulant method (if expanded in Taylor
374 // series, one recovers NewG() in new new cumulant method)
375
376 TVector2 normQ = NormQ(pFlowSelect) ;
377
378 return exp(2*Zx*normQ.X() + 2*Zy*normQ.Y());
379}
380//-----------------------------------------------------------------------
381Double_t AliFlowEvent::SumWeightSquare(AliFlowSelection* pFlowSelect) const
382{
383 // Return sum of weights^2 for the specified Selection (used for normalization)
384
385 if(fDone) { return fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
386
387 int selN = pFlowSelect->Sel();
388 int harN = pFlowSelect->Har();
389 Double_t sumOfWeightSqr = 0;
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);
398 sumOfWeightSqr += phiWgt*phiWgt;
399 }
400 }
401 if(sumOfWeightSqr < 0.) { return Mult(pFlowSelect) ; }
402
403 return sumOfWeightSqr;
404}
405//-------------------------------------------------------------
406Double_t AliFlowEvent::WgtMultQ4(AliFlowSelection* pFlowSelect) const
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//-------------------------------------------------------------
442Double_t AliFlowEvent::WgtMultQ6(AliFlowSelection* pFlowSelect) const
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) *
506 for(int harN=0;harN<AliFlowConstants::kHars;harN++) { pFlowTrack->SetSelect(harN,0) ; }
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());
528 float pt = pFlowTrack->Pt();
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
535 for(int selN = 1; selN < AliFlowConstants::kSels; selN++) // not even consider the 0th selection (no cut applied there)
536 {
537 // min. TPC hits required
538 if(pFlowSelect->NhitsCut(selN) && (pFlowTrack->FitPtsTPC()<pFlowSelect->NhitsCut(selN))) continue ;
539
540 for(int harN = 0; harN < AliFlowConstants::kHars; harN++)
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
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 ;
546
547 pFlowTrack->SetSelect(harN, selN) ; // if cuts defined (low<high) && track is in the range -> Set [har][sel] Flag ON
548
549 if(AliFlowConstants::fgDebug)
550 {
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 ;
561 }
562 }
563 }
564 }
565}
566//-------------------------------------------------------------
567void AliFlowEvent::SetPids()
568{
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}.
571
572 const Int_t kCode[] = {11,13,211,321,2212,10010020} ;
573 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
574 {
575 AliFlowTrack* pFlowTrack ;
576 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
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)
584 for(Int_t nP=0;nP<AliFlowConstants::kPid;nP++)
585 {
586 if(bayPid[nP]>pidMax) { maxN = nP ; pidMax = bayPid[nP] ; }
587 }
588 Int_t pdgCode = TMath::Sign(kCode[maxN],pFlowTrack->Charge()) ;
589 pFlowTrack->SetMostLikelihoodPID(pdgCode) ;
590 }
591}
592//-------------------------------------------------------------
593void AliFlowEvent::MakeSubEvents() const
594{
595 // Make random or eta sub-events
596
597 if(fgEtaSubs) { MakeEtaSubEvents() ; }
598 else { MakeRndSubEvents() ; }
599}
600//-------------------------------------------------------------
601void AliFlowEvent::MakeRndSubEvents() const
602{
603 // Make random subevents
604
605 int eventMult[AliFlowConstants::kHars][AliFlowConstants::kSels] = {{0}};
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) ;
613 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
614 {
615 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
616 {
617 if(pFlowTrack->Select(harN, selN)) { eventMult[harN][selN]++ ; }
618 }
619 }
620 }
621 // loop to set the SubEvent member
622 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
623 {
624 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
625 {
626 int subEventMult = eventMult[harN][selN] / AliFlowConstants::kSubs;
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//-------------------------------------------------------------
648void AliFlowEvent::MakeEtaSubEvents() const
649{
650 // Make subevents for positive and negative eta
651 // (when done, fgEtaSubs flag setted to kTRUE).
652
653 int harN, selN = 0;
654 // loop to set the SubEvent member
655 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
656 {
657 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
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
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 ;
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//-----------------------------------------------------------------------
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 {
731 cent = AliFlowConstants::fgCentNorm ;
732 //if centrality classes are not defined, does it now (with CentNorm & MaxMult)
733 if(cent[AliFlowConstants::kCents-1] <= 1)
734 {
735 for(Int_t ic=0;ic<AliFlowConstants::kCents;ic++)
736 {
737 cent[ic] *= AliFlowConstants::fgMaxMult ;
738 if(AliFlowConstants::fgDebug) { cout << "Centrality[" << ic << "] = " << cent[ic] << " . " << endl ; }
739 }
740 }
741 }
742 else if((RunID() != -1) && (CenterOfMassEnergy() == 5500.))
743 {
744 cent = (Float_t*)AliFlowConstants::fgCent0 ;
745 }
746 else // other definition of centrality are possible ...
747 {
748 cent = (Float_t*)AliFlowConstants::fgCent0 ;
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
761 if(AliFlowConstants::fgDebug) { cout << " * Centrality Class : " << fCentrality << " . " << endl ; }
762}
763//-----------------------------------------------------------------------
764TVector AliFlowEvent::Bayesian() const
765{
766 // Returns bayesian array of particle abundances (from AliFlowConstants::)
767
768 TVector bayes(AliFlowConstants::kPid) ;
769 for(Int_t i=0;i<AliFlowConstants::kPid;i++) { bayes[i] = AliFlowConstants::fgBayesian[i] ; }
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?)
782 cout << "# PtWgt Saturation = " << AliFlowConstants::fgPtWgtSaturation << endl;
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//-----------------------------------------------------------------------
799Int_t AliFlowEvent::MultEta() const
800{
801 // Returns the multiplicity in the interval |eta|<(AliFlowConstants::fgEetaMid), used
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) ;
809 if((pFlowTrack->Charge()) && (TMath::Abs(pFlowTrack->Eta())<AliFlowConstants::fgEtaMid)) { goodtracks++ ; }
810 }
811 return goodtracks ;
812}
813//-----------------------------------------------------------------------
814Int_t AliFlowEvent::UncorrNegMult(Float_t eta) const
815{
816 // Negative multiplicity in the interval (-eta..eta)
817 // (default is AliFlowConstants::fgEetaGood = 0.9)
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)
833 // (default is AliFlowConstants::fgEetaGood = 0.9)
834
835 Int_t posMult = 0 ;
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//-----------------------------------------------------------------------
846TVector3 AliFlowEvent::VertexPos() const
847{
848 // Returns primary vertex position as a TVector3
849
850 Float_t vertex[3] ;
851 VertexPos(vertex) ;
852 return TVector3(vertex) ;
853}
854//-----------------------------------------------------------------------
855void AliFlowEvent::MakeAll()
856{
857 // calculates all quantities in 1 shoot
858
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] ;
863 // -
864 int selN, harN, subN ;
865 for(selN=0;selN<AliFlowConstants::kSels;selN++)
866 {
867 for(harN=0;harN<AliFlowConstants::kHars;harN++)
868 {
869 mQx[selN][harN] = 0. ;
870 mQy[selN][harN] = 0. ;
871 fMult[selN][harN] = 0 ;
872 fSumOfWeightSqr[selN][harN] = 0. ;
873 for(subN=0;subN<AliFlowConstants::kSubs;subN++)
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();
892 for(selN=0;selN<AliFlowConstants::kSels;selN++)
893 {
894 for(harN=0;harN<AliFlowConstants::kHars;harN++)
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]++ ;
904 for(subN=0;subN<AliFlowConstants::kSubs;subN++)
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
918 for(selN=0;selN<AliFlowConstants::kSels;selN++)
919 {
920 for(harN=0;harN<AliFlowConstants::kHars;harN++)
921 {
922 fQ[selN][harN].Set(mQx[selN][harN],mQy[selN][harN]) ;
923 for(subN=0;subN<AliFlowConstants::kSubs;subN++) { fQSub[subN][selN][harN].Set(mQxSub[subN][selN][harN],mQySub[subN][selN][harN]) ; }
924 }
925 }
926
927 fDone = kTRUE ;
928}
929//-----------------------------------------------------------------------