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