PWG2/FLOW/Attic/AliFlowEvent.cxx

   1 //////////////////////////////////////////////////////////////////////
   2 //
   3 // $Id: AliFlowEvent.cxx 18618 2007-05-16 15:38:22Z snelling $
   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>
  52 using namespace std; //required for resolving the 'cout' symbol
  53
  54 // - Flags & Sets
  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)
  62
  63 ClassImp(AliFlowEvent)
  64 //-----------------------------------------------------------
  65 AliFlowEvent::AliFlowEvent(Int_t lenght):
  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)
  75 {
  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
  79
  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.; }
  82
  83  // Make a new track collection
  84  fTrackCollection =  new TClonesArray("AliFlowTrack",lenght) ; fTrackCollection->BypassStreamer(kTRUE) ;
  85  fV0Collection    =  new TClonesArray("AliFlowV0",lenght)    ; fV0Collection->BypassStreamer(kTRUE)    ;
  86
  87  // Set Weights Arrays to 1 (default)
  88  for(int nS=0;nS<AliFlowConstants::kSels;nS++)
  89  {
  90   for(int nH=0;nH<AliFlowConstants::kHars;nH++)
  91   {
  92    for(int nP=0;nP<AliFlowConstants::kPhiBins;nP++)
  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  }
 103  //for(int nH=0;nH<AliFlowConstants::kHars;nH++) { fExtPsi[nH] = 0. ; fExtRes[nH] = 0. ; }
 104
 105  // The Expected particles abundance is taken directly from AliFlowConstants::fgBayesian[] (see Bayesian P.Id.)
 106
 107 }
 108 //-----------------------------------------------------------
 109 AliFlowEvent::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 //-----------------------------------------------------------
 140 AliFlowEvent::~AliFlowEvent()
 141 {
 142  // Default distructor: deletes the ObjArrays.
 143
 144  if(fTrackCollection) { delete fTrackCollection ; } // fTrackCollection->Delete() ;
 145  if(fV0Collection)    { delete fV0Collection ; }    // fV0Collection->Delete()    ;
 146 }
 147 //-------------------------------------------------------------
 148
 149 //-------------------------------------------------------------
 150 Double_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() ;
 156  int n = (int)((phi/(2*TMath::Pi()))*AliFlowConstants::kPhiBins);
 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
 172  return phiWgt ;
 173 }
 174 //-------------------------------------------------------------
 175 Double_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
 179  if(selN>AliFlowConstants::kSels) { selN = 0 ; }
 180  bool oddHar = (harN+1) % 2 ;
 181  Double_t phiWgt = 1. ;
 182  if(PtWgt())
 183  {
 184   Double_t pt = (Double_t)pFlowTrack->Pt();
 185   if(pt < AliFlowConstants::fgPtWgtSaturation) { phiWgt *= pt ; }
 186   else                            { phiWgt *= AliFlowConstants::fgPtWgtSaturation ; } // pt weighting going constant
 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 //-------------------------------------------------------------
 196 Double_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
 199  // (it simply rerurns PhiWeightRaw() * Weight()). If fgNoWgt = kTRUE, returns +/-1 ,
 200  // basing on Sign(eta), for odd harmonics .
 201
 202  if(fgNoWgt) // no weights (but +/- according to eta)
 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);
 210  if(AliFlowConstants::fgDebug) { cout << "[PhiWeight]: phiWgtRaw = " << phiWgtRaw << " , weight = " << weight << " , eta = " << pFlowTrack->Eta() << endl ; }
 211
 212  return phiWgtRaw * weight;
 213 }
 214 //-------------------------------------------------------------
 215 Int_t AliFlowEvent::Mult(AliFlowSelection* pFlowSelect) const
 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 //-------------------------------------------------------------
 237 Float_t AliFlowEvent::MeanPt(AliFlowSelection* pFlowSelect) const
 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 //-------------------------------------------------------------
 260 TVector2 AliFlowEvent::Q(AliFlowSelection* pFlowSelect) const
 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) ;
 288    if(AliFlowConstants::fgDebug) { cout << itr << " phi = " << phi << " ,  wgt = " << phiWgt << endl ; }
 289   }
 290  }
 291  mQ.Set(mQx, mQy);
 292
 293  return mQ;
 294 }
 295 //-------------------------------------------------------------
 296 Float_t AliFlowEvent::Psi(AliFlowSelection* pFlowSelect) const
 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 //-------------------------------------------------------------
 313 TVector2 AliFlowEvent::NormQ(AliFlowSelection* pFlowSelect) const
 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()] ;
 320   double sumOfWeightSqr = fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ;
 321   if(sumOfWeightSqr) { mQ /= TMath::Sqrt(sumOfWeightSqr) ; }
 322   else { mQ.Set(0.,0.) ; }
 323   return mQ ;
 324  }
 325  // -
 326  TVector2 mQ ;
 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) ;
 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);
 340    sumOfWeightSqr += phiWgt*phiWgt;
 341
 342    float phi = pFlowTrack->Phi();
 343    mQx += phiWgt * cos(phi * order);
 344    mQy += phiWgt * sin(phi * order);
 345   }
 346  }
 347  if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
 348  else { mQ.Set(0.,0.); }
 349
 350  return mQ;
 351 }
 352 //-------------------------------------------------------------
 353 Float_t AliFlowEvent::OldQ(AliFlowSelection* pFlowSelect) const
 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) ;
 362  double sumOfWeightSqr = 0 ;
 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
 372    sumOfWeightSqr += phiWgt*phiWgt;
 373    float phi = pFlowTrack->Phi();
 374    mQx += phiWgt * cos(phi * order);
 375    mQy += phiWgt * sin(phi * order);
 376   }
 377  }
 378  if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
 379  else { mQ.Set(0.,0.); }
 380
 381  return mQ.Mod();
 382 }
 383 //-----------------------------------------------------------------------
 384 Double_t AliFlowEvent::NewG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
 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 //-----------------------------------------------------------------------
 410 Double_t AliFlowEvent::OldG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
 411 {
 412  // Generating function for the old cumulant method (if expanded in Taylor
 413  // series, one recovers NewG() in new new cumulant method)
 414
 415  TVector2 normQ = NormQ(pFlowSelect) ;
 416
 417  return exp(2*Zx*normQ.X() + 2*Zy*normQ.Y());
 418 }
 419 //-----------------------------------------------------------------------
 420 Double_t AliFlowEvent::SumWeightSquare(AliFlowSelection* pFlowSelect) const
 421 {
 422  // Return sum of weights^2 for the specified Selection (used for normalization)
 423
 424  if(fDone) { return fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
 425
 426  int selN = pFlowSelect->Sel();
 427  int harN = pFlowSelect->Har();
 428  Double_t sumOfWeightSqr = 0;
 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);
 437    sumOfWeightSqr += phiWgt*phiWgt;
 438   }
 439  }
 440  if(sumOfWeightSqr < 0.) { return Mult(pFlowSelect) ; }
 441
 442  return sumOfWeightSqr;
 443 }
 444 //-------------------------------------------------------------
 445 Double_t AliFlowEvent::WgtMultQ4(AliFlowSelection* pFlowSelect) const
 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 //-------------------------------------------------------------
 481 Double_t AliFlowEvent::WgtMultQ6(AliFlowSelection* pFlowSelect) const
 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 //-------------------------------------------------------------
 520 void 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) *
 545   for(int harN=0;harN<AliFlowConstants::kHars;harN++) { pFlowTrack->SetSelect(harN,0) ; }
 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());
 567   float  pt  = pFlowTrack->Pt();
 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
 574   for(int selN = 1; selN < AliFlowConstants::kSels; selN++)               // not even consider the 0th selection (no cut applied there)
 575   {
 576   // min. TPC hits required
 577    if(pFlowSelect->NhitsCut(selN) && (pFlowTrack->FitPtsTPC()<pFlowSelect->NhitsCut(selN))) continue ;
 578
 579    for(int harN = 0; harN < AliFlowConstants::kHars; harN++)
 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
 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 ;
 585
 586     pFlowTrack->SetSelect(harN, selN) ;  // if cuts defined (low<high) && track is in the range -> Set [har][sel] Flag ON
 587
 588     if(AliFlowConstants::fgDebug)
 589     {
 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 ;
 600     }
 601    }
 602   }
 603  }
 604 }
 605 //-------------------------------------------------------------
 606 void AliFlowEvent::SetPids()
 607 {
 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}.
 610
 611  const Int_t kCode[] = {11,13,211,321,2212,10010020} ;
 612  for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++)
 613  {
 614   AliFlowTrack* pFlowTrack ;
 615   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
 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)
 623   for(Int_t nP=0;nP<AliFlowConstants::kPid;nP++)
 624   {
 625    if(bayPid[nP]>pidMax) { maxN = nP ; pidMax = bayPid[nP] ; }
 626   }
 627   Int_t pdgCode = TMath::Sign(kCode[maxN],pFlowTrack->Charge()) ;
 628   pFlowTrack->SetMostLikelihoodPID(pdgCode) ;
 629  }
 630 }
 631 //-------------------------------------------------------------
 632 void AliFlowEvent::MakeSubEvents() const
 633 {
 634  // Make random or eta sub-events
 635
 636  if(fgEtaSubs == 1)        { MakeEtaSubEvents() ; }
 637  else if(fgEtaSubs == -1)  { MakeChrSubEvents() ; }
 638  else if(fgEtaSubs == 0)   { MakeRndSubEvents() ; }
 639 }
 640 //-------------------------------------------------------------
 641 void AliFlowEvent::MakeRndSubEvents() const
 642 {
 643  // Make random subevents
 644
 645  int eventMult[AliFlowConstants::kHars][AliFlowConstants::kSels] = {{0}};
 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) ;
 653   for (selN = 0; selN < AliFlowConstants::kSels; selN++)
 654   {
 655    for (harN = 0; harN < AliFlowConstants::kHars; harN++)
 656    {
 657     if(pFlowTrack->Select(harN, selN)) { eventMult[harN][selN]++ ; }
 658    }
 659   }
 660  }
 661  // loop to set the SubEvent member
 662  for (selN = 0; selN < AliFlowConstants::kSels; selN++)
 663  {
 664   for (harN = 0; harN < AliFlowConstants::kHars; harN++)
 665   {
 666    int subEventMult = eventMult[harN][selN] / AliFlowConstants::kSubs;
 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 //-------------------------------------------------------------
 688 void AliFlowEvent::MakeEtaSubEvents() const
 689 {
 690  // Make subevents for positive and negative eta
 691
 692  int harN, selN = 0;
 693  // loop to set the SubEvent member
 694  for (selN = 0; selN < AliFlowConstants::kSels; selN++)
 695  {
 696   for (harN = 0; harN < AliFlowConstants::kHars; harN++)
 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 //-------------------------------------------------------------
 713 void 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 //-------------------------------------------------------------
 738 void 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
 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 ;
 775 }
 776 //-----------------------------------------------------------------------
 777 UInt_t AliFlowEvent::Centrality()
 778 {
 779  // Returns the Centrality Class as stored
 780
 781  if(fCentrality < 0)  { SetCentrality() ; }
 782  return fCentrality ;
 783 }
 784 //-----------------------------------------------------------------------
 785 void 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  {
 795   cent = AliFlowConstants::fgCentNorm ;
 796   //if centrality classes are not defined, does it now (with CentNorm & MaxMult)
 797   if(cent[AliFlowConstants::kCents-1] <= 1)
 798   {
 799    for(Int_t ic=0;ic<AliFlowConstants::kCents;ic++)
 800    {
 801     cent[ic] *= AliFlowConstants::fgMaxMult ;
 802     if(AliFlowConstants::fgDebug) { cout << "Centrality[" << ic << "] = " << cent[ic] << " . " << endl ; }
 803    }
 804   }
 805  }
 806  else if((RunID() != -1) && (CenterOfMassEnergy() == 5500.))
 807  {
 808   cent = (Float_t*)AliFlowConstants::fgCent0 ;
 809  }
 810  else // other definition of centrality are possible ...
 811  {
 812   cent = (Float_t*)AliFlowConstants::fgCent0 ;
 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
 825  if(AliFlowConstants::fgDebug) { cout << " * Centrality Class :  " << fCentrality << " . " << endl ; }
 826 }
 827 //-----------------------------------------------------------------------
 828 TVector AliFlowEvent::Bayesian() const
 829 {
 830  // Returns bayesian array of particle abundances (from AliFlowConstants::)
 831
 832  TVector bayes(AliFlowConstants::kPid) ;
 833  for(Int_t i=0;i<AliFlowConstants::kPid;i++) { bayes[i] = AliFlowConstants::fgBayesian[i] ; }
 834  return bayes ;
 835 }
 836 //-----------------------------------------------------------------------
 837 void 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?)
 846   cout << "#    PtWgt Saturation =  " << AliFlowConstants::fgPtWgtSaturation << endl;
 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 //-----------------------------------------------------------------------
 863 Int_t AliFlowEvent::MultEta() const
 864 {
 865  // Returns the multiplicity in the interval |eta|<(AliFlowConstants::fgEetaMid), used
 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) ;
 873   if((pFlowTrack->Charge()) && (TMath::Abs(pFlowTrack->Eta())<AliFlowConstants::fgEtaMid)) { goodtracks++ ; }
 874  }
 875  return goodtracks ;
 876 }
 877 //-----------------------------------------------------------------------
 878 Int_t AliFlowEvent::UncorrNegMult(Float_t eta)  const
 879 {
 880  // Negative multiplicity in the interval (-eta..eta)
 881  // (default is  AliFlowConstants::fgEetaGood = 0.9)
 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++ ; }
 889   //delete pFlowTrack ;
 890  }
 891  return negMult;
 892 }
 893 //-----------------------------------------------------------------------
 894 Int_t AliFlowEvent::UncorrPosMult(Float_t eta)  const
 895 {
 896  // Positive multiplicity in the interval (-eta..eta)
 897  // (default is  AliFlowConstants::fgEetaGood = 0.9)
 898
 899  Int_t posMult = 0 ;
 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++ ; }
 905   //delete pFlowTrack ;
 906  }
 907  return posMult;
 908 }
 909 //-----------------------------------------------------------------------
 910 TVector3 AliFlowEvent::VertexPos() const
 911 {
 912  // Returns primary vertex position as a TVector3
 913
 914  Float_t vertex[3] ;
 915  VertexPos(vertex) ;
 916  return TVector3(vertex) ;
 917 }
 918 //-----------------------------------------------------------------------
 919 void AliFlowEvent::MakeAll()
 920 {
 921  // calculates all quantities in 1 shoot
 922
 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] ;
 927  // -
 928  int selN, harN, subN ;
 929  for(selN=0;selN<AliFlowConstants::kSels;selN++)
 930  {
 931   for(harN=0;harN<AliFlowConstants::kHars;harN++)
 932   {
 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++)
 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();
 956   for(selN=0;selN<AliFlowConstants::kSels;selN++)
 957   {
 958    for(harN=0;harN<AliFlowConstants::kHars;harN++)
 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]++ ;
 968      for(subN=0;subN<AliFlowConstants::kSubs;subN++)
 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
 982  for(selN=0;selN<AliFlowConstants::kSels;selN++)
 983  {
 984   for(harN=0;harN<AliFlowConstants::kHars;harN++)
 985   {
 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]) ; }
 988   }
 989  }
 990
 991  fDone = kTRUE ;
 992 }
 993 //-----------------------------------------------------------------------