PWG2/FLOW/AliFlowEvent.cxx

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