1 //////////////////////////////////////////////////////////////////////
5 // Author: Emanuele Simili
7 //////////////////////////////////////////////////////////////////////
9 //_____________________________________________________________
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).
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
38 #include "AliFlowEvent.h"
39 #include "AliFlowTrack.h"
40 #include "AliFlowV0.h"
41 #include "AliFlowSelection.h"
42 #include "AliFlowConstants.h"
47 #include "TClonesArray.h"
52 using namespace std; //required for resolving the 'cout' symbol
55 Bool_t AliFlowEvent::fgPtWgt = kFALSE ; // gives pT-based weights
56 Bool_t AliFlowEvent::fgEtaWgt = kFALSE ; // gives eta-based weights
57 Bool_t AliFlowEvent::fgOnePhiWgt = kTRUE ; // kTRUE --> ENABLEs SINGLE WEIGHT ISTOGRAM , kFALSE --> ENABLEs 3 WEIGHT ISTOGRAMS
58 Bool_t AliFlowEvent::fgNoWgt = kFALSE ; // No Weight is used
59 Bool_t AliFlowEvent::fgCustomRespFunc = kFALSE ; // custom "detector response function" is used for P.Id (see AliFlowTrack)
60 // - Eta Sub-Events (later used to calculate the resolution)
61 Int_t AliFlowEvent::fgEtaSubs = 0 ; // makes random subevents (0 = random , 1 = eta , -1 = charged)
63 ClassImp(AliFlowEvent)
64 //-----------------------------------------------------------
65 AliFlowEvent::AliFlowEvent(Int_t lenght):
72 fTrackCollection(0x0),
76 // not-Default constructor: initializes the ObjArray of FlowTracks and FlowV0s,
77 // cleans the internal variables, sets all the weights to 1, sets default flags.
78 // USE THIS WHEN CREATING ALIFLOWEVENTS
80 for(int zz=0;zz<3;zz++) { fZDCenergy[zz] = 0. ; }
81 for(int i=0;i<AliFlowConstants::kHars;i++) { fExtPsi[i] = 0. ; fExtRes[i] = 0.; }
83 // Make a new track collection
84 fTrackCollection = new TClonesArray("AliFlowTrack",lenght) ; fTrackCollection->BypassStreamer(kTRUE) ;
85 fV0Collection = new TClonesArray("AliFlowV0",lenght) ; fV0Collection->BypassStreamer(kTRUE) ;
87 // Set Weights Arrays to 1 (default)
88 for(int nS=0;nS<AliFlowConstants::kSels;nS++)
90 for(int nH=0;nH<AliFlowConstants::kHars;nH++)
92 for(int nP=0;nP<AliFlowConstants::kPhiBins;nP++)
94 // enable this with: SetOnePhiWgt()
95 fPhiWgt[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgt[nS][nH][nP] << endl ;
96 // enable these with: SetFirstLastPhiWgt()
97 fPhiWgtPlus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtPlus[nS][nH][nP] << endl ;
98 fPhiWgtMinus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtMinus[nS][nH][nP] << endl ;
99 fPhiWgtCross[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtCross[nS][nH][nP] << endl ;
103 //for(int nH=0;nH<AliFlowConstants::kHars;nH++) { fExtPsi[nH] = 0. ; fExtRes[nH] = 0. ; }
105 // The Expected particles abundance is taken directly from AliFlowConstants::fgBayesian[] (see Bayesian P.Id.)
108 //-----------------------------------------------------------
109 AliFlowEvent::AliFlowEvent():
116 fTrackCollection(0x0),
120 // Default constructor : USE THIS WHEN READING ALIFLOWEVENTS
122 // Set Weights Arrays to 1 (default)
123 for(int nS=0;nS<AliFlowConstants::kSels;nS++)
125 for(int nH=0;nH<AliFlowConstants::kHars;nH++)
127 for(int nP=0;nP<AliFlowConstants::kPhiBins;nP++)
129 // enable this with: SetOnePhiWgt()
130 fPhiWgt[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgt[nS][nH][nP] << endl ;
131 // enable these with: SetFirstLastPhiWgt()
132 fPhiWgtPlus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtPlus[nS][nH][nP] << endl ;
133 fPhiWgtMinus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtMinus[nS][nH][nP] << endl ;
134 fPhiWgtCross[nS][nH][nP] = 1. ; // cout << nS << nH << nP << " val " << fPhiWgtCross[nS][nH][nP] << endl ;
139 //-----------------------------------------------------------
140 AliFlowEvent::~AliFlowEvent()
142 // Default distructor: deletes the ObjArrays.
144 if(fTrackCollection) { delete fTrackCollection ; } // fTrackCollection->Delete() ;
145 if(fV0Collection) { delete fV0Collection ; } // fV0Collection->Delete() ;
147 //-------------------------------------------------------------
149 //-------------------------------------------------------------
150 Double_t AliFlowEvent::PhiWeightRaw(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
152 // Weight for making the event plane isotropic in the lab.
154 float phi = pFlowTrack->Phi() ; if(phi < 0.) { phi += 2*TMath::Pi() ; }
155 Double_t eta = (Double_t)pFlowTrack->Eta() ;
156 int n = (int)((phi/(2*TMath::Pi()))*AliFlowConstants::kPhiBins);
158 Double_t phiWgt = 1. ;
161 phiWgt = (Double_t)fPhiWgt[selN][harN][n]; //cout << "Selection " << selN << " ; Harmonic " << harN << " ; PhiBin " << n << " - Wgt = " << phiWgt << endl ;
163 else if(FirstLastPhiWgt())
165 float zFirstPoint = pFlowTrack->ZFirstPoint(); // float zLastPoint = pFlowTrack->ZLastPoint();
167 if (zFirstPoint > 0. && eta > 0.) { phiWgt = (Double_t)fPhiWgtPlus[selN][harN][n] ; }
168 else if(zFirstPoint < 0. && eta < 0.) { phiWgt = (Double_t)fPhiWgtMinus[selN][harN][n] ; }
169 else { phiWgt = (Double_t)fPhiWgtCross[selN][harN][n] ; }
174 //-------------------------------------------------------------
175 Double_t AliFlowEvent::Weight(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
177 // Weight for enhancing the resolution (eta gives sign +/- for Odd Harmonics)
179 if(selN>AliFlowConstants::kSels) { selN = 0 ; }
180 bool oddHar = (harN+1) % 2 ;
181 Double_t phiWgt = 1. ;
184 Double_t pt = (Double_t)pFlowTrack->Pt();
185 if(pt < AliFlowConstants::fgPtWgtSaturation) { phiWgt *= pt ; }
186 else { phiWgt *= AliFlowConstants::fgPtWgtSaturation ; } // pt weighting going constant
188 Double_t eta = (Double_t)pFlowTrack->Eta();
189 Double_t etaAbs = TMath::Abs(eta);
190 if(EtaWgt() && oddHar) { phiWgt *= etaAbs ; }
191 if(oddHar && eta < 0.) { phiWgt *= -1. ; }
195 //-------------------------------------------------------------
196 Double_t AliFlowEvent::PhiWeight(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
198 // Weight for making the event plane isotropic in the lab and enhancing the resolution
199 // (it simply rerurns PhiWeightRaw() * Weight()). If fgNoWgt = kTRUE, returns +/-1 ,
200 // basing on Sign(eta), for odd harmonics .
202 if(fgNoWgt) // no weights (but +/- according to eta)
204 bool oddHar = (harN+1) % 2 ;
205 if(oddHar) { return TMath::Sign((Double_t)1.,(Double_t)pFlowTrack->Eta()) ; }
206 else { return (Double_t)1. ; }
208 Double_t phiWgtRaw = PhiWeightRaw(selN, harN, pFlowTrack);
209 Double_t weight = Weight(selN, harN, pFlowTrack);
210 if(AliFlowConstants::fgDebug) { cout << "[PhiWeight]: phiWgtRaw = " << phiWgtRaw << " , weight = " << weight << " , eta = " << pFlowTrack->Eta() << endl ; }
212 return phiWgtRaw * weight;
214 //-------------------------------------------------------------
215 Int_t AliFlowEvent::Mult(AliFlowSelection* pFlowSelect) const
217 // Multiplicity of tracks in the specified Selection
221 int sub = pFlowSelect->Sub() ;
222 if(sub<0) { return fMult[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
223 else { return fMultSub[sub][pFlowSelect->Sel()][pFlowSelect->Har()] ; }
228 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
230 AliFlowTrack* pFlowTrack ;
231 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
232 if(pFlowSelect->Select(pFlowTrack)) { mult++ ; }
236 //-------------------------------------------------------------
237 Float_t AliFlowEvent::MeanPt(AliFlowSelection* pFlowSelect) const
239 // Mean pt of tracks in the specified Selection
241 Double_t meanPt = 0. ;
245 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
247 AliFlowTrack* pFlowTrack ;
248 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
249 if (pFlowSelect->Select(pFlowTrack))
251 sumPt += pFlowTrack->Pt();
255 if(mult) { meanPt = sumPt/(float)mult ; }
259 //-------------------------------------------------------------
260 TVector2 AliFlowEvent::Q(AliFlowSelection* pFlowSelect) const
262 // Event plane vector for the specified Selection
266 int sub = pFlowSelect->Sub() ;
267 if(sub<0) { return fQ[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
268 else { return fQSub[sub][pFlowSelect->Sel()][pFlowSelect->Har()] ; }
272 Double_t mQx=0. , mQy=0. ;
273 int selN = pFlowSelect->Sel() ;
274 int harN = pFlowSelect->Har() ;
275 double order = (double)(harN + 1) ;
278 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
280 AliFlowTrack* pFlowTrack ;
281 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
282 if(pFlowSelect->Select(pFlowTrack))
284 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
285 float phi = pFlowTrack->Phi();
286 mQx += phiWgt * cos(phi * order) ;
287 mQy += phiWgt * sin(phi * order) ;
288 if(AliFlowConstants::fgDebug) { cout << itr << " phi = " << phi << " , wgt = " << phiWgt << endl ; }
295 //-------------------------------------------------------------
296 Float_t AliFlowEvent::Psi(AliFlowSelection* pFlowSelect) const
298 // Event plane angle for the specified Selection
300 int harN = pFlowSelect->Har() ;
301 float order = (float)(harN + 1) ;
304 TVector2 mQ = Q(pFlowSelect);
305 if(mQ.Mod()) // if vector is not 0
307 psi= mQ.Phi() / order ;
308 if (psi < 0.) { psi += 2*TMath::Pi() / order ; }
312 //-------------------------------------------------------------
313 TVector2 AliFlowEvent::NormQ(AliFlowSelection* pFlowSelect) const
315 // Normalized Q = Q/sqrt(sum of weights^2) for the specified Selection
319 TVector2 mQ = fQ[pFlowSelect->Sel()][pFlowSelect->Har()] ;
320 double sumOfWeightSqr = fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ;
321 if(sumOfWeightSqr) { mQ /= TMath::Sqrt(sumOfWeightSqr) ; }
322 else { mQ.Set(0.,0.) ; }
327 Double_t mQx=0. , mQy=0. ;
328 double sumOfWeightSqr = 0 ;
329 int selN = pFlowSelect->Sel() ;
330 int harN = pFlowSelect->Har() ;
331 double order = (double)(harN + 1) ;
333 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
335 AliFlowTrack* pFlowTrack ;
336 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
337 if (pFlowSelect->Select(pFlowTrack))
339 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
340 sumOfWeightSqr += phiWgt*phiWgt;
342 float phi = pFlowTrack->Phi();
343 mQx += phiWgt * cos(phi * order);
344 mQy += phiWgt * sin(phi * order);
347 if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
348 else { mQ.Set(0.,0.); }
352 //-------------------------------------------------------------
353 Float_t AliFlowEvent::OldQ(AliFlowSelection* pFlowSelect) const
355 // Magnitude of normalized Q vector (without pt or eta weighting) for the specified Selection
358 Double_t mQx = 0. , mQy = 0. ;
359 int selN = pFlowSelect->Sel() ;
360 int harN = pFlowSelect->Har() ;
361 double order = (double)(harN + 1) ;
362 double sumOfWeightSqr = 0 ;
365 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
367 AliFlowTrack* pFlowTrack ;
368 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
369 if(pFlowSelect->Select(pFlowTrack))
371 double phiWgt = PhiWeightRaw(selN, harN, pFlowTrack); // Raw
372 sumOfWeightSqr += phiWgt*phiWgt;
373 float phi = pFlowTrack->Phi();
374 mQx += phiWgt * cos(phi * order);
375 mQy += phiWgt * sin(phi * order);
378 if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
379 else { mQ.Set(0.,0.); }
383 //-----------------------------------------------------------------------
384 Double_t AliFlowEvent::NewG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
386 // Generating function for the new cumulant method. Eq. 3 in the Practical Guide
388 int selN = pFlowSelect->Sel();
389 int harN = pFlowSelect->Har();
390 double order = (double)(harN + 1);
392 double mult = (double)Mult(pFlowSelect);
396 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
398 AliFlowTrack* pFlowTrack ;
399 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
400 if (pFlowSelect->Select(pFlowTrack))
402 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
403 float phi = pFlowTrack->Phi();
404 theG *= (1. + (phiWgt/mult) * (2.* Zx * cos(phi * order) + 2.* Zy * sin(phi * order) ) );
409 //-----------------------------------------------------------------------
410 Double_t AliFlowEvent::OldG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
412 // Generating function for the old cumulant method (if expanded in Taylor
413 // series, one recovers NewG() in new new cumulant method)
415 TVector2 normQ = NormQ(pFlowSelect) ;
417 return exp(2*Zx*normQ.X() + 2*Zy*normQ.Y());
419 //-----------------------------------------------------------------------
420 Double_t AliFlowEvent::SumWeightSquare(AliFlowSelection* pFlowSelect) const
422 // Return sum of weights^2 for the specified Selection (used for normalization)
424 if(fDone) { return fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
426 int selN = pFlowSelect->Sel();
427 int harN = pFlowSelect->Har();
428 Double_t sumOfWeightSqr = 0;
430 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
432 AliFlowTrack* pFlowTrack ;
433 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
434 if (pFlowSelect->Select(pFlowTrack))
436 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
437 sumOfWeightSqr += phiWgt*phiWgt;
440 if(sumOfWeightSqr < 0.) { return Mult(pFlowSelect) ; }
442 return sumOfWeightSqr;
444 //-------------------------------------------------------------
445 Double_t AliFlowEvent::WgtMultQ4(AliFlowSelection* pFlowSelect) const
447 // Used only for the old cumulant method, for getting q4 when weight is on.
448 // Replace multiplicity in Eq.(74b) by this quantity when weight is on.
449 // This is derived based on (A4) in the old cumulant paper.
451 int selN = pFlowSelect->Sel();
452 int harN = pFlowSelect->Har();
454 double theMeanWj4 = 0.;
455 double theMeanWj2 = 0.;
456 double theSumOfWgtSqr = 0;
460 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
462 AliFlowTrack* pFlowTrack ;
463 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
464 if (pFlowSelect->Select(pFlowTrack))
466 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
467 phiWgtSq = phiWgt*phiWgt;
468 theSumOfWgtSqr += phiWgtSq;
469 theMeanWj4 += phiWgtSq*phiWgtSq;
473 if (theMult <= 0.) return theMult;
475 theMeanWj4 /= theMult;
476 theMeanWj2 = theSumOfWgtSqr / theMult;
478 return (theSumOfWgtSqr*theSumOfWgtSqr)/(theMult*(-theMeanWj4+2*theMeanWj2*theMeanWj2));
480 //-------------------------------------------------------------
481 Double_t AliFlowEvent::WgtMultQ6(AliFlowSelection* pFlowSelect) const
483 // Used only for the old cumulant method. For getting q6 when weight is on.
484 // Replace multiplicity in Eq.(74c) by this quantity when weight is on.
485 // This is derived based on (A4) in the old cumulant paper.
487 int selN = pFlowSelect->Sel();
488 int harN = pFlowSelect->Har();
490 double theMeanWj6 = 0.;
491 double theMeanWj4 = 0.;
492 double theMeanWj2 = 0.;
493 double theSumOfWgtSqr = 0;
497 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
499 AliFlowTrack* pFlowTrack ;
500 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
501 if (pFlowSelect->Select(pFlowTrack))
503 double phiWgt = PhiWeight(selN, harN, pFlowTrack);
504 phiWgtSq = phiWgt*phiWgt;
505 theSumOfWgtSqr += phiWgtSq;
506 theMeanWj4 += phiWgtSq*phiWgtSq;
507 theMeanWj6 += phiWgtSq*phiWgtSq*phiWgtSq;
511 if (theMult <= 0.) return theMult*theMult;
513 theMeanWj6 /= theMult;
514 theMeanWj4 /= theMult;
515 theMeanWj2 = theSumOfWgtSqr / theMult;
517 return 4.*(theSumOfWgtSqr*theSumOfWgtSqr*theSumOfWgtSqr)/(theMult*(theMeanWj6-9.*theMeanWj2*theMeanWj4+12.*theMeanWj2*theMeanWj2*theMeanWj2));
519 //-------------------------------------------------------------
520 void AliFlowEvent::SetSelections(AliFlowSelection* pFlowSelect)
522 // Sets the selection of tracks used in the Reaction Plane calculation
523 // for the specific Harmonic and Selection - this does not cut trow away
524 // tracks from the event, neither exclude them from flow study. See also
525 // the AliFlowSelection class.
526 // Strategy of Selection:
527 // For the specific Harmonic and Selection, IF cuts
528 // are defined (such that low<high) and IF the track satisfies them, THEN
529 // the respective track's flag (bool AliFlowTrack::mSelection[har][sel])
530 // is set kTRUE so that the track, from now on, will be included in the
531 // R.P. determination for that selection.
532 // If NO cuts are defined -> ALL Flags are setted kTRUE (all particle
533 // used for all the Reaction Planes -> no difference in Psi[har][sel]).
534 // -------------------------------------------------------------------
535 // The first selection (all harmonics) is set kTRUE : no conditions.
538 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
540 AliFlowTrack* pFlowTrack ;
541 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
542 pFlowTrack->ResetSelection() ; // this re-sets all the mSelection flags to 0
544 // * this sets all the selection n.[0] flag kTRUE (all harmonics) *
545 for(int harN=0;harN<AliFlowConstants::kHars;harN++) { pFlowTrack->SetSelect(harN,0) ; }
547 // Track need to be Constrainable
548 if(pFlowSelect->ConstrainCut() && !pFlowTrack->IsConstrainable()) continue ;
550 // PID - gets rid of the track with AliFlowTrack::Pid() != AliFlowSelection::Pid() (if there)
551 if(pFlowSelect->Pid()[0] != '\0')
553 if(strstr(pFlowSelect->Pid(), "h")!=0)
555 int charge = pFlowTrack->Charge() ;
556 if(strcmp("h+",pFlowSelect->Pid())==0 && charge != 1) continue;
557 if(strcmp("h-",pFlowSelect->Pid())==0 && charge != -1) continue;
562 strcpy(pid, pFlowTrack->Pid());
563 if(strstr(pid, pFlowSelect->Pid())==0) continue;
566 double eta = (double)(pFlowTrack->Eta());
567 float pt = pFlowTrack->Pt();
568 float gDca = pFlowTrack->Dca() ;
570 // Global DCA - gets rid of the track with DCA outside the range
571 if((pFlowSelect->DcaGlobalCutHi()>pFlowSelect->DcaGlobalCutLo()) && (gDca<pFlowSelect->DcaGlobalCutLo() || gDca>pFlowSelect->DcaGlobalCutHi())) continue ;
573 // Pt & Eta - this is done differently for different Harmonic & Selection
574 for(int selN = 1; selN < AliFlowConstants::kSels; selN++) // not even consider the 0th selection (no cut applied there)
576 // min. TPC hits required
577 if(pFlowSelect->NhitsCut(selN) && (pFlowTrack->FitPtsTPC()<pFlowSelect->NhitsCut(selN))) continue ;
579 for(int harN = 0; harN < AliFlowConstants::kHars; harN++)
581 // Eta - gets rid of the track with Eta outside the range
582 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 ;
583 // Pt - gets rid of the track with Pt outside the range
584 if((pFlowSelect->PtCutHi(harN%2,selN)>pFlowSelect->PtCutLo(harN%2,selN)) && (pt<pFlowSelect->PtCutLo(harN%2,selN) || pt>pFlowSelect->PtCutHi(harN%2,selN))) continue ;
586 pFlowTrack->SetSelect(harN, selN) ; // if cuts defined (low<high) && track is in the range -> Set [har][sel] Flag ON
588 if(AliFlowConstants::fgDebug)
592 cout << " n. " << itr << " , pFlowTrack->PrintSelection() = " ; pFlowTrack->PrintSelection() ;
593 if(pFlowSelect->Pid()[0] != '\0') { cout << " track: pid " << pFlowTrack->Pid() << " = "<< pFlowSelect->Pid() << endl ; }
594 cout << " track: dca " << pFlowSelect->DcaGlobalCutLo() << " < " << gDca << " < " << pFlowSelect->DcaGlobalCutHi() << endl ;
595 cout << " track: eta " << pFlowSelect->EtaCutLo(harN,selN) << " < |" << eta << "| < " << pFlowSelect->EtaCutHi(harN,selN) << endl ;
596 cout << " track: pT " << pFlowSelect->PtCutLo(harN,selN) << " < " << pt << " < " << pFlowSelect->PtCutHi(harN,selN) << endl ;
597 cout << " pFlowTrack->PrintSelection() = " ; pFlowTrack->PrintSelection() ;
599 // cout << " selN " << selN << " , harN " << harN%2 << " - si" << endl ;
605 //-------------------------------------------------------------
606 void AliFlowEvent::SetPids()
608 // Re-sets the tracks P.id. (using the current AliFlowConstants::fgBayesian[] array).
609 // To use the Raw P.Id (just the detector response function), set fgBayesian[] = {1,1,1,1,1,1}.
611 const Int_t kCode[] = {11,13,211,321,2212,10010020} ;
612 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
614 AliFlowTrack* pFlowTrack ;
615 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
617 Float_t bayPid[AliFlowConstants::kPid] ;
618 if(fgCustomRespFunc) { pFlowTrack->PidProbsC(bayPid) ; }
619 else { pFlowTrack->PidProbs(bayPid) ; }
621 Int_t maxN = 2 ; // if No id. -> then is a Pi
622 Float_t pidMax = bayPid[2] ; // (if all equal, Pi probability was the first)
623 for(Int_t nP=0;nP<AliFlowConstants::kPid;nP++)
625 if(bayPid[nP]>pidMax) { maxN = nP ; pidMax = bayPid[nP] ; }
627 Int_t pdgCode = TMath::Sign(kCode[maxN],pFlowTrack->Charge()) ;
628 pFlowTrack->SetMostLikelihoodPID(pdgCode) ;
631 //-------------------------------------------------------------
632 void AliFlowEvent::MakeSubEvents() const
634 // Make random or eta sub-events
636 if(fgEtaSubs == 1) { MakeEtaSubEvents() ; }
637 else if(fgEtaSubs == -1) { MakeChrSubEvents() ; }
638 else if(fgEtaSubs == 0) { MakeRndSubEvents() ; }
640 //-------------------------------------------------------------
641 void AliFlowEvent::MakeRndSubEvents() const
643 // Make random subevents
645 int eventMult[AliFlowConstants::kHars][AliFlowConstants::kSels] = {{0}};
646 int harN, selN, subN = 0;
648 // loop to count the total number of tracks for each selection
649 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
651 AliFlowTrack* pFlowTrack ;
652 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
653 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
655 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
657 if(pFlowTrack->Select(harN, selN)) { eventMult[harN][selN]++ ; }
661 // loop to set the SubEvent member
662 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
664 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
666 int subEventMult = eventMult[harN][selN] / AliFlowConstants::kSubs;
671 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
673 AliFlowTrack* pFlowTrack ;
674 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
675 if(pFlowTrack->Select(harN, selN))
677 pFlowTrack->SetSubevent(harN, selN, subN);
679 if (countN % subEventMult == 0.) { subN++ ; }
687 //-------------------------------------------------------------
688 void AliFlowEvent::MakeEtaSubEvents() const
690 // Make subevents for positive and negative eta
693 // loop to set the SubEvent member
694 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
696 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
698 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
700 AliFlowTrack* pFlowTrack ;
701 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
702 if(pFlowTrack->Select(harN, selN))
704 float eta = pFlowTrack->Eta();
705 if (eta > 0.) { pFlowTrack->SetSubevent(harN, selN, 0) ; }
706 else { pFlowTrack->SetSubevent(harN, selN, 1) ; }
712 //-------------------------------------------------------------
713 void AliFlowEvent::MakeChrSubEvents() const
715 // Make subevents for positive and negative charged tracks
718 // loop to set the SubEvent member
719 for (selN = 0; selN < AliFlowConstants::kSels; selN++)
721 for (harN = 0; harN < AliFlowConstants::kHars; harN++)
723 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
725 AliFlowTrack* pFlowTrack ;
726 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
727 if(pFlowTrack->Select(harN, selN))
729 float charge = pFlowTrack->Charge();
730 if (charge > 0.) { pFlowTrack->SetSubevent(harN, selN, 0) ; }
731 else { pFlowTrack->SetSubevent(harN, selN, 1) ; }
737 //-------------------------------------------------------------
738 void AliFlowEvent::RandomShuffle()
740 // Randomly re-shuffles the tracks in the array; if a track is not
741 // primary, the reference carried by the reconstructed mother (in
742 // the v0 array) is changed accordingly.
744 return ; // ...at the moment is disabled ... //
747 // UInt_t imax = fTrackCollection->GetEntries() ;
748 // TRandom* rnd = new TRandom(0) ; // SetSeed(0) ;
749 // TClonesArray* newTrackCollection = new TClonesArray("AliFlowTrack",imax) ;
751 // // re-arranges the ObjArray (TrackCollection())
752 // for(UInt_t itr=0;itr<imax;itr++)
754 // AliFlowTrack* pFlowTrack ;
755 // pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
757 // UInt_t rndNumber = rnd->Integer(imax) ;
758 // Bool_t put = kFALSE ;
761 // if(!newTrackCollection->AddrAt(rndNumber))
763 // ... new(newTrackCollection[rndNumber]) AliFlowTrack(*pFlowTrack) ;
764 // put = kTRUE ; tot++ ;
765 // if(AliFlowConstants::fgDebug) { cout << " " << itr << " --> " << rndNumber << endl ; }
769 // rndNumber++ ; if(rndNumber>=imax) { rndNumber -= imax ; }
773 // if(AliFlowConstants::fgDebug) { cout << "* RandomShuffle() : " << tot << "/" << imax << " flow tracks have been shuffled " << endl ; }
774 // fTrackCollection = newTrackCollection ;
776 //-----------------------------------------------------------------------
777 UInt_t AliFlowEvent::Centrality()
779 // Returns the Centrality Class as stored
781 if(fCentrality < 0) { SetCentrality() ; }
784 //-----------------------------------------------------------------------
785 void AliFlowEvent::SetCentrality()
787 // Sets the Centrality Classes basing on Multiplicity at mid rapidity,
788 // ... (ZDC information can be added) .
791 Int_t tracks = MultEta() ;
795 cent = AliFlowConstants::fgCentNorm ;
796 //if centrality classes are not defined, does it now (with CentNorm & MaxMult)
797 if(cent[AliFlowConstants::kCents-1] <= 1)
799 for(Int_t ic=0;ic<AliFlowConstants::kCents;ic++)
801 cent[ic] *= AliFlowConstants::fgMaxMult ;
802 if(AliFlowConstants::fgDebug) { cout << "Centrality[" << ic << "] = " << cent[ic] << " . " << endl ; }
806 else if((RunID() != -1) && (CenterOfMassEnergy() == 5500.))
808 cent = (Float_t*)AliFlowConstants::fgCent0 ;
810 else // other definition of centrality are possible ...
812 cent = (Float_t*)AliFlowConstants::fgCent0 ;
814 if (tracks < cent[0]) { fCentrality = 0; }
815 else if (tracks < cent[1]) { fCentrality = 1; }
816 else if (tracks < cent[2]) { fCentrality = 2; }
817 else if (tracks < cent[3]) { fCentrality = 3; }
818 else if (tracks < cent[4]) { fCentrality = 4; }
819 else if (tracks < cent[5]) { fCentrality = 5; }
820 else if (tracks < cent[6]) { fCentrality = 6; }
821 else if (tracks < cent[7]) { fCentrality = 7; }
822 else if (tracks < cent[8]) { fCentrality = 8; }
823 else { fCentrality = 9; }
825 if(AliFlowConstants::fgDebug) { cout << " * Centrality Class : " << fCentrality << " . " << endl ; }
827 //-----------------------------------------------------------------------
828 TVector AliFlowEvent::Bayesian() const
830 // Returns bayesian array of particle abundances (from AliFlowConstants::)
832 TVector bayes(AliFlowConstants::kPid) ;
833 for(Int_t i=0;i<AliFlowConstants::kPid;i++) { bayes[i] = AliFlowConstants::fgBayesian[i] ; }
836 //-----------------------------------------------------------------------
837 void AliFlowEvent::PrintFlagList() const
839 // Prints the list of selection cuts ( called in AliFlowInterface::Finish() )
841 cout << "#######################################################" << endl;
842 cout << "# Weighting and Striping:" << endl;
845 cout << "# PtWgt = kTRUE " << endl ; // (also for output of PhiWgt file?)
846 cout << "# PtWgt Saturation = " << AliFlowConstants::fgPtWgtSaturation << endl;
850 cout << "# PtWgt = kFALSE" << endl;
854 cout << "# EtaWgt = kTRUE " << endl ; // (also for output of PhiWgt file for odd harmonics?)
858 cout << "# EtaWgt = kFALSE" << endl;
860 cout << "#######################################################" << endl;
862 //-----------------------------------------------------------------------
863 Int_t AliFlowEvent::MultEta() const
865 // Returns the multiplicity in the interval |eta|<(AliFlowConstants::fgEetaMid), used
866 // for centrality measurement (see centrality classes in fCentrality) .
868 Int_t goodtracks = 0 ;
869 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
871 AliFlowTrack* pFlowTrack ;
872 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
873 if((pFlowTrack->Charge()) && (TMath::Abs(pFlowTrack->Eta())<AliFlowConstants::fgEtaMid)) { goodtracks++ ; }
877 //-----------------------------------------------------------------------
878 Int_t AliFlowEvent::UncorrNegMult(Float_t eta) const
880 // Negative multiplicity in the interval (-eta..eta)
881 // (default is AliFlowConstants::fgEetaGood = 0.9)
884 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
886 AliFlowTrack* pFlowTrack ;
887 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
888 if((pFlowTrack->Charge()<0) && (TMath::Abs(pFlowTrack->Eta())<TMath::Abs(eta))) { negMult++ ; }
889 //delete pFlowTrack ;
893 //-----------------------------------------------------------------------
894 Int_t AliFlowEvent::UncorrPosMult(Float_t eta) const
896 // Positive multiplicity in the interval (-eta..eta)
897 // (default is AliFlowConstants::fgEetaGood = 0.9)
900 for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
902 AliFlowTrack* pFlowTrack ;
903 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
904 if((pFlowTrack->Charge()>0) && (TMath::Abs(pFlowTrack->Eta())<TMath::Abs(eta))) { posMult++ ; }
905 //delete pFlowTrack ;
909 //-----------------------------------------------------------------------
910 TVector3 AliFlowEvent::VertexPos() const
912 // Returns primary vertex position as a TVector3
916 return TVector3(vertex) ;
918 //-----------------------------------------------------------------------
919 void AliFlowEvent::MakeAll()
921 // calculates all quantities in 1 shoot
923 Double_t mQx[AliFlowConstants::kSels][AliFlowConstants::kHars] ;
924 Double_t mQy[AliFlowConstants::kSels][AliFlowConstants::kHars] ;
925 Double_t mQxSub[AliFlowConstants::kSubs][AliFlowConstants::kSels][AliFlowConstants::kHars] ;
926 Double_t mQySub[AliFlowConstants::kSubs][AliFlowConstants::kSels][AliFlowConstants::kHars] ;
928 int selN, harN, subN ;
929 for(selN=0;selN<AliFlowConstants::kSels;selN++)
931 for(harN=0;harN<AliFlowConstants::kHars;harN++)
933 mQx[selN][harN] = 0. ;
934 mQy[selN][harN] = 0. ;
935 fMult[selN][harN] = 0 ;
936 fSumOfWeightSqr[selN][harN] = 0. ;
937 for(subN=0;subN<AliFlowConstants::kSubs;subN++)
939 mQxSub[subN][selN][harN] = 0. ;
940 mQySub[subN][selN][harN] = 0. ;
941 fMultSub[subN][selN][harN] = 0 ;
951 for(itr=0;itr<TrackCollection()->GetEntries();itr++)
953 AliFlowTrack* pFlowTrack ;
954 pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
955 phi = pFlowTrack->Phi();
956 for(selN=0;selN<AliFlowConstants::kSels;selN++)
958 for(harN=0;harN<AliFlowConstants::kHars;harN++)
960 order = (double)(harN+1) ;
961 if(pFlowTrack->Select(harN,selN))
963 phiWgt = PhiWeight(selN,harN,pFlowTrack) ;
964 fSumOfWeightSqr[selN][harN] += phiWgt*phiWgt ;
965 mQx[selN][harN] += phiWgt * cos(phi * order) ;
966 mQy[selN][harN] += phiWgt * sin(phi * order) ;
967 fMult[selN][harN]++ ;
968 for(subN=0;subN<AliFlowConstants::kSubs;subN++)
970 if(pFlowTrack->Select(harN,selN,subN))
972 mQxSub[subN][selN][harN] += phiWgt * cos(phi * order) ;
973 mQySub[subN][selN][harN] += phiWgt * sin(phi * order) ;
974 fMultSub[subN][selN][harN]++ ;
982 for(selN=0;selN<AliFlowConstants::kSels;selN++)
984 for(harN=0;harN<AliFlowConstants::kHars;harN++)
986 fQ[selN][harN].Set(mQx[selN][harN],mQy[selN][harN]) ;
987 for(subN=0;subN<AliFlowConstants::kSubs;subN++) { fQSub[subN][selN][harN].Set(mQxSub[subN][selN][harN],mQySub[subN][selN][harN]) ; }
993 //-----------------------------------------------------------------------