PWG/FLOW/Base/AliFlowEventSimple.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*****************************************************************
  17   AliFlowEventSimple: A simple event
  18   for flow analysis
  19
  20   origin: Naomi van der Kolk (kolk@nikhef.nl)
  21           Ante Bilandzic     (anteb@nikhef.nl)
  22           Raimond Snellings  (Raimond.Snellings@nikhef.nl)
  23   mods:   Mikolaj Krzewicki  (mikolaj.krzewicki@cern.ch)
  24 *****************************************************************/
  25
  26 #include "Riostream.h"
  27 #include "TObjArray.h"
  28 #include "TFile.h"
  29 #include "TList.h"
  30 #include "TTree.h"
  31 #include "TParticle.h"
  32 #include "TMath.h"
  33 #include "TH1F.h"
  34 #include "TH1D.h"
  35 #include "TF1.h"
  36 #include "TProfile.h"
  37 #include "TParameter.h"
  38 #include "TBrowser.h"
  39 #include "AliFlowVector.h"
  40 #include "AliFlowTrackSimple.h"
  41 #include "AliFlowTrackSimpleCuts.h"
  42 #include "AliFlowEventSimple.h"
  43 #include "TRandom.h"
  44
  45 ClassImp(AliFlowEventSimple)
  46
  47 //-----------------------------------------------------------------------
  48 AliFlowEventSimple::AliFlowEventSimple():
  49   fTrackCollection(NULL),
  50   fReferenceMultiplicity(0),
  51   fNumberOfTracks(0),
  52   fNumberOfRPs(0),
  53   fMCReactionPlaneAngle(0.),
  54   fMCReactionPlaneAngleIsSet(kFALSE),
  55   fAfterBurnerPrecision(0.001),
  56   fUserModified(kFALSE),
  57   fNumberOfTracksWrap(NULL),
  58   fNumberOfRPsWrap(NULL),
  59   fMCReactionPlaneAngleWrap(NULL)
  60 {
  61   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
  62 }
  63
  64 //-----------------------------------------------------------------------
  65 AliFlowEventSimple::AliFlowEventSimple( Int_t n,
  66                                         ConstructionMethod method,
  67                                         TF1* ptDist,
  68                                         Double_t phiMin,
  69                                         Double_t phiMax,
  70                                         Double_t etaMin,
  71                                         Double_t etaMax):
  72   fTrackCollection(new TObjArray(n)),
  73   fReferenceMultiplicity(0),
  74   fNumberOfTracks(0),
  75   fNumberOfRPs(0),
  76   fMCReactionPlaneAngle(0.),
  77   fMCReactionPlaneAngleIsSet(kFALSE),
  78   fAfterBurnerPrecision(0.001),
  79   fUserModified(kFALSE),
  80   fNumberOfTracksWrap(NULL),
  81   fNumberOfRPsWrap(NULL),
  82   fMCReactionPlaneAngleWrap(NULL)
  83 {
  84   //ctor
  85   // if second argument is set to AliFlowEventSimple::kGenerate
  86   // it generates n random tracks with given Pt distribution
  87   // (a sane default is provided), phi and eta are uniform
  88
  89   if (method==kGenerate)
  90     Generate(n,ptDist,phiMin,phiMax,etaMin,etaMax);
  91 }
  92
  93 //-----------------------------------------------------------------------
  94 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
  95   TObject(anEvent),
  96   fTrackCollection((TObjArray*)(anEvent.fTrackCollection)->Clone()),
  97   fReferenceMultiplicity(anEvent.fReferenceMultiplicity),
  98   fNumberOfTracks(anEvent.fNumberOfTracks),
  99   fNumberOfRPs(anEvent.fNumberOfRPs),
 100   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
 101   fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
 102   fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
 103   fUserModified(anEvent.fUserModified),
 104   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
 105   fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
 106   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 107 {
 108   //copy constructor
 109 }
 110
 111 //-----------------------------------------------------------------------
 112 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
 113 {
 114   //assignment operator
 115   if (&anEvent==this) return *this; //check self-assignment
 116   if (fTrackCollection) fTrackCollection->Delete();
 117   delete fTrackCollection;
 118   fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
 119   fReferenceMultiplicity = anEvent.fReferenceMultiplicity;
 120   fNumberOfTracks = anEvent.fNumberOfTracks;
 121   fNumberOfRPs = anEvent.fNumberOfRPs;
 122   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
 123   fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
 124   fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
 125   fUserModified=anEvent.fUserModified;
 126   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
 127   fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
 128   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
 129   return *this;
 130 }
 131
 132 //-----------------------------------------------------------------------
 133 AliFlowEventSimple::~AliFlowEventSimple()
 134 {
 135   //destructor
 136   if (fTrackCollection) fTrackCollection->Delete();
 137   delete fTrackCollection;
 138   delete fNumberOfTracksWrap;
 139   delete fNumberOfRPsWrap;
 140   delete fMCReactionPlaneAngleWrap;
 141 }
 142
 143 //-----------------------------------------------------------------------
 144 void AliFlowEventSimple::Generate(Int_t nParticles,
 145                                   TF1* ptDist,
 146                                   Double_t phiMin,
 147                                   Double_t phiMax,
 148                                   Double_t etaMin,
 149                                   Double_t etaMax)
 150 {
 151   //generate nParticles random tracks uniform in phi and eta
 152   //according to the specified pt distribution
 153   if (!ptDist)
 154   {
 155     static TF1 ptdistribution("ptSpectra","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
 156     ptDist=&ptdistribution;
 157   }
 158
 159   for (Int_t i=0; i<nParticles; i++)
 160   {
 161     AliFlowTrackSimple* track = new AliFlowTrackSimple();
 162     track->SetPhi( gRandom->Uniform(phiMin,phiMax) );
 163     track->SetEta( gRandom->Uniform(etaMin,etaMax) );
 164     track->SetPt( ptDist->GetRandom() );
 165     track->SetCharge( (gRandom->Uniform()-0.5<0)?-1:1 );
 166     AddTrack(track);
 167   }
 168   fMCReactionPlaneAngle=gRandom->Uniform(0.0,TMath::TwoPi());
 169   fMCReactionPlaneAngleIsSet=kTRUE;
 170   SetUserModified();
 171 }
 172
 173 //-----------------------------------------------------------------------
 174 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 175 {
 176   //get track i from collection
 177   if (i>=fNumberOfTracks) return NULL;
 178   AliFlowTrackSimple* pTrack = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i)) ;
 179   return pTrack;
 180 }
 181
 182 //-----------------------------------------------------------------------
 183 void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
 184 {
 185   //add a track, delete the old one if necessary
 186   if (fNumberOfTracks < fTrackCollection->GetEntriesFast())
 187   {
 188     TObject* o = fTrackCollection->At(fNumberOfTracks);
 189     if (o) delete o;
 190   }
 191   fTrackCollection->AddAtAndExpand(track,fNumberOfTracks);
 192   fNumberOfTracks++;
 193 }
 194
 195 //-----------------------------------------------------------------------
 196 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 197 {
 198   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 199   Double_t dQX = 0.;
 200   Double_t dQY = 0.;
 201   AliFlowVector vQ;
 202   vQ.Set(0.,0.);
 203
 204   Int_t iOrder = n;
 205   Double_t sumOfWeights = 0.;
 206   Double_t dPhi = 0.;
 207   Double_t dPt = 0.;
 208   Double_t dEta = 0.;
 209   Double_t dWeight = 1.;
 210
 211   AliFlowTrackSimple* pTrack = NULL;
 212
 213   Int_t nBinsPhi = 0;
 214   Double_t dBinWidthPt = 0.;
 215   Double_t dPtMin = 0.;
 216   Double_t dBinWidthEta = 0.;
 217   Double_t dEtaMin = 0.;
 218
 219   Double_t wPhi = 1.; // weight Phi
 220   Double_t wPt = 1.;  // weight Pt
 221   Double_t wEta = 1.; // weight Eta
 222
 223   TH1F *phiWeights = NULL;
 224   TH1D *ptWeights  = NULL;
 225   TH1D *etaWeights = NULL;
 226
 227   if(weightsList)
 228   {
 229     if(usePhiWeights)
 230     {
 231       phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
 232       if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
 233     }
 234     if(usePtWeights)
 235     {
 236       ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 237       if(ptWeights)
 238       {
 239         dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 240         dPtMin = (ptWeights->GetXaxis())->GetXmin();
 241       }
 242     }
 243     if(useEtaWeights)
 244     {
 245       etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 246       if(etaWeights)
 247       {
 248         dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 249         dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 250       }
 251     }
 252   } // end of if(weightsList)
 253
 254   // loop over tracks
 255   for(Int_t i=0; i<fNumberOfTracks; i++)
 256   {
 257     pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
 258     if(pTrack)
 259     {
 260       if(pTrack->InRPSelection())
 261       {
 262         dPhi = pTrack->Phi();
 263         dPt  = pTrack->Pt();
 264         dEta = pTrack->Eta();
 265               dWeight = pTrack->Weight();
 266
 267         // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 268         if(phiWeights && nBinsPhi)
 269         {
 270           wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
 271         }
 272         // determine v'(pt) weight:
 273         if(ptWeights && dBinWidthPt)
 274         {
 275           wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 276         }
 277         // determine v'(eta) weight:
 278         if(etaWeights && dBinWidthEta)
 279         {
 280           wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 281         }
 282
 283         // building up the weighted Q-vector:
 284         dQX += dWeight*wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
 285         dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
 286
 287         // weighted multiplicity:
 288         sumOfWeights += dWeight*wPhi*wPt*wEta;
 289
 290       } // end of if (pTrack->InRPSelection())
 291     } // end of if (pTrack)
 292     else
 293     {
 294       cerr << "no particle!!!"<<endl;
 295     }
 296   } // loop over particles
 297
 298   vQ.Set(dQX,dQY);
 299   vQ.SetMult(sumOfWeights);
 300
 301   return vQ;
 302
 303 }
 304
 305 //-----------------------------------------------------------------------
 306 void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 307 {
 308
 309   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 310   Double_t dQX = 0.;
 311   Double_t dQY = 0.;
 312
 313   Int_t iOrder = n;
 314   Double_t sumOfWeights = 0.;
 315   Double_t dPhi = 0.;
 316   Double_t dPt  = 0.;
 317   Double_t dEta = 0.;
 318   Double_t dWeight = 1.;
 319
 320   AliFlowTrackSimple* pTrack = NULL;
 321
 322   Int_t    iNbinsPhiSub0 = 0;
 323   Int_t    iNbinsPhiSub1 = 0;
 324   Double_t dBinWidthPt = 0.;
 325   Double_t dPtMin      = 0.;
 326   Double_t dBinWidthEta= 0.;
 327   Double_t dEtaMin     = 0.;
 328
 329   Double_t dWphi = 1.;  // weight Phi
 330   Double_t dWpt  = 1.;  // weight Pt
 331   Double_t dWeta = 1.;  // weight Eta
 332
 333   TH1F* phiWeightsSub0 = NULL;
 334   TH1F* phiWeightsSub1 = NULL;
 335   TH1D* ptWeights  = NULL;
 336   TH1D* etaWeights = NULL;
 337
 338   if(weightsList)
 339   {
 340     if(usePhiWeights)
 341     {
 342       phiWeightsSub0 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub0"));
 343       if(phiWeightsSub0) {
 344         iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
 345       }
 346       phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
 347       if(phiWeightsSub1) {
 348         iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
 349       }
 350     }
 351     if(usePtWeights)
 352     {
 353       ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 354       if(ptWeights)
 355       {
 356         dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 357         dPtMin = (ptWeights->GetXaxis())->GetXmin();
 358       }
 359     }
 360     if(useEtaWeights)
 361     {
 362       etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 363       if(etaWeights)
 364       {
 365         dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 366         dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 367       }
 368     }
 369   } // end of if(weightsList)
 370
 371   //loop over the two subevents
 372   for (Int_t s=0; s<2; s++)
 373   {
 374     // loop over tracks
 375     for(Int_t i=0; i<fNumberOfTracks; i++)
 376     {
 377       pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
 378       if(pTrack)
 379       {
 380         if(pTrack->InRPSelection())
 381         {
 382           if (pTrack->InSubevent(s))
 383           {
 384             dPhi    = pTrack->Phi();
 385             dPt     = pTrack->Pt();
 386             dEta    = pTrack->Eta();
 387             dWeight = pTrack->Weight();
 388
 389             // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 390             //subevent 0
 391             if(s == 0)  {
 392               if(phiWeightsSub0 && iNbinsPhiSub0)  {
 393                 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi()));
 394                 //use the phi value at the center of the bin
 395                 dPhi  = phiWeightsSub0->GetBinCenter(phiBin);
 396                 dWphi = phiWeightsSub0->GetBinContent(phiBin);
 397               }
 398             }
 399             //subevent 1
 400             else if (s == 1) {
 401               if(phiWeightsSub1 && iNbinsPhiSub1) {
 402                 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi()));
 403                 //use the phi value at the center of the bin
 404                 dPhi  = phiWeightsSub1->GetBinCenter(phiBin);
 405                 dWphi = phiWeightsSub1->GetBinContent(phiBin);
 406               }
 407             }
 408
 409             // determine v'(pt) weight:
 410             if(ptWeights && dBinWidthPt)
 411             {
 412               dWpt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 413             }
 414
 415             // determine v'(eta) weight:
 416             if(etaWeights && dBinWidthEta)
 417             {
 418               dWeta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 419             }
 420
 421             // building up the weighted Q-vector:
 422             dQX += dWeight*dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
 423             dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
 424
 425             // weighted multiplicity:
 426             sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
 427
 428           } // end of subevent
 429         } // end of if (pTrack->InRPSelection())
 430       } // end of if (pTrack)
 431       else
 432       {
 433         cerr << "no particle!!!"<<endl;
 434       }
 435     } // loop over particles
 436     Qarray[s].Set(dQX,dQY);
 437     Qarray[s].SetMult(sumOfWeights);
 438     //reset
 439     sumOfWeights = 0.;
 440     dQX = 0.;
 441     dQY = 0.;
 442   }
 443
 444 }
 445
 446
 447 //-----------------------------------------------------------------------
 448 void AliFlowEventSimple::Print(Option_t *option) const
 449 {
 450   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
 451   //             ===============================================
 452   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
 453   printf( "Class.Print Name = %s, #tracks= %d, Number of RPs= %d, MC EventPlaneAngle= %f\n",
 454           GetName(),fNumberOfTracks, fNumberOfRPs, fMCReactionPlaneAngle );
 455
 456   TString optionstr(option);
 457   if (!optionstr.Contains("all")) return;
 458   if (fTrackCollection)
 459   {
 460     fTrackCollection->Print(option);
 461   }
 462   else
 463   {
 464     printf( "Empty track collection \n");
 465   }
 466 }
 467
 468 //-----------------------------------------------------------------------
 469 void AliFlowEventSimple::Browse(TBrowser *b)
 470 {
 471   //browse in TBrowser
 472   if (!b) return;
 473   if (!fNumberOfTracksWrap)
 474   {
 475     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
 476     b->Add(fNumberOfTracksWrap);
 477   }
 478   if (!fNumberOfRPsWrap)
 479   {
 480     fNumberOfRPsWrap = new TParameter<int>("fNumberOfRPs", fNumberOfRPs);
 481     b->Add(fNumberOfRPsWrap);
 482   }
 483   if (!fMCReactionPlaneAngleWrap)
 484   {
 485     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
 486     b->Add( fMCReactionPlaneAngleWrap);
 487   }
 488   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
 489 }
 490
 491 //-----------------------------------------------------------------------
 492 AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
 493                                         const AliFlowTrackSimpleCuts* rpCuts,
 494                                         const AliFlowTrackSimpleCuts* poiCuts):
 495   fTrackCollection(NULL),
 496   fReferenceMultiplicity(0),
 497   fNumberOfTracks(0),
 498   fNumberOfRPs(0),
 499   fMCReactionPlaneAngle(0.),
 500   fMCReactionPlaneAngleIsSet(kFALSE),
 501   fAfterBurnerPrecision(0.001),
 502   fUserModified(kFALSE),
 503   fNumberOfTracksWrap(NULL),
 504   fNumberOfRPsWrap(NULL),
 505   fMCReactionPlaneAngleWrap(NULL)
 506 {
 507   //constructor, fills the event from a TTree of kinematic.root files
 508   //applies RP and POI cuts, tags the tracks
 509
 510   Int_t numberOfInputTracks = inputTree->GetEntries() ;
 511   fTrackCollection = new TObjArray(numberOfInputTracks/2);
 512
 513   TParticle* pParticle = new TParticle();
 514   inputTree->SetBranchAddress("Particles",&pParticle);
 515
 516   Int_t iSelParticlesPOI = 0;
 517
 518   for (Int_t i=0; i<numberOfInputTracks; i++)
 519   {
 520     inputTree->GetEntry(i);   //get input particle
 521
 522     if (!pParticle) continue; //no particle
 523     if (!pParticle->IsPrimary()) continue;
 524
 525     Bool_t rpOK = rpCuts->PassesCuts(pParticle);
 526     Bool_t poiOK = poiCuts->PassesCuts(pParticle);
 527
 528     if (rpOK || poiOK)
 529     {
 530       AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);
 531
 532       //marking the particles used for int. flow:
 533       if(rpOK)
 534       {
 535         pTrack->SetForRPSelection(kTRUE);
 536         fNumberOfRPs++;
 537       }
 538       //marking the particles used for diff. flow:
 539       if(poiOK)
 540       {
 541         pTrack->SetForPOISelection(kTRUE);
 542         iSelParticlesPOI++;
 543       }
 544       //adding a particles which were used either for int. or diff. flow to the list
 545       AddTrack(pTrack);
 546     }
 547   }//for i
 548   delete pParticle;
 549 }
 550
 551 //_____________________________________________________________________________
 552 void AliFlowEventSimple::CloneTracks(Int_t n)
 553 {
 554   //clone every track n times to add non-flow
 555   if (n<=0) return; //no use to clone stuff zero or less times
 556   Int_t ntracks = fNumberOfTracks;
 557   fTrackCollection->Expand((n+1)*fNumberOfTracks);
 558   for (Int_t i=0; i<n; i++)
 559   {
 560     for (Int_t itrack=0; itrack<ntracks; itrack++)
 561     {
 562       AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
 563       if (!track) continue;
 564       AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
 565     }
 566   }
 567   SetUserModified();
 568 }
 569
 570 //_____________________________________________________________________________
 571 void AliFlowEventSimple::ResolutionPt(Double_t res)
 572 {
 573   //smear pt of all tracks by gaussian with sigma=res
 574   for (Int_t i=0; i<fNumberOfTracks; i++)
 575   {
 576     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 577     if (track) track->ResolutionPt(res);
 578   }
 579   SetUserModified();
 580 }
 581
 582 //_____________________________________________________________________________
 583 void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
 584                                             Double_t etaMaxA,
 585                                             Double_t etaMinB,
 586                                             Double_t etaMaxB )
 587 {
 588   //Flag two subevents in given eta ranges
 589   for (Int_t i=0; i<fNumberOfTracks; i++)
 590   {
 591     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 592     if (!track) continue;
 593     track->ResetSubEventTags();
 594     Double_t eta=track->Eta();
 595     if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
 596     if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
 597   }
 598 }
 599
 600 //_____________________________________________________________________________
 601 void AliFlowEventSimple::TagSubeventsByCharge()
 602 {
 603   //Flag two subevents in given eta ranges
 604   for (Int_t i=0; i<fNumberOfTracks; i++)
 605   {
 606     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 607     if (!track) continue;
 608     track->ResetSubEventTags();
 609     Int_t charge=track->Charge();
 610     if (charge<0) track->SetForSubevent(0);
 611     if (charge>0) track->SetForSubevent(1);
 612   }
 613 }
 614
 615 //_____________________________________________________________________________
 616 void AliFlowEventSimple::AddV1( Double_t v1 )
 617 {
 618   //add v2 to all tracks wrt the reaction plane angle
 619   for (Int_t i=0; i<fNumberOfTracks; i++)
 620   {
 621     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 622     if (track) track->AddV1(v1, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 623   }
 624   SetUserModified();
 625 }
 626
 627 //_____________________________________________________________________________
 628 void AliFlowEventSimple::AddV2( Double_t v2 )
 629 {
 630   //add v2 to all tracks wrt the reaction plane angle
 631   for (Int_t i=0; i<fNumberOfTracks; i++)
 632   {
 633     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 634     if (track) track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 635   }
 636   SetUserModified();
 637 }
 638
 639 //_____________________________________________________________________________
 640 void AliFlowEventSimple::AddV3( Double_t v3 )
 641 {
 642   //add v3 to all tracks wrt the reaction plane angle
 643   for (Int_t i=0; i<fNumberOfTracks; i++)
 644   {
 645     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 646     if (track) track->AddV3(v3, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 647   }
 648   SetUserModified();
 649 }
 650
 651 //_____________________________________________________________________________
 652 void AliFlowEventSimple::AddV4( Double_t v4 )
 653 {
 654   //add v4 to all tracks wrt the reaction plane angle
 655   for (Int_t i=0; i<fNumberOfTracks; i++)
 656   {
 657     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 658     if (track) track->AddV4(v4, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 659   }
 660   SetUserModified();
 661 }
 662
 663 //_____________________________________________________________________________
 664 void AliFlowEventSimple::AddV5( Double_t v5 )
 665 {
 666   //add v4 to all tracks wrt the reaction plane angle
 667   for (Int_t i=0; i<fNumberOfTracks; i++)
 668   {
 669     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 670     if (track) track->AddV5(v5, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 671   }
 672   SetUserModified();
 673 }
 674
 675 //_____________________________________________________________________________
 676 void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5,
 677                                   Double_t rp1, Double_t rp2, Double_t rp3, Double_t rp4, Double_t rp5 )
 678 {
 679   //add flow to all tracks wrt the reaction plane angle, for all harmonic separate angle
 680   for (Int_t i=0; i<fNumberOfTracks; i++)
 681   {
 682     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 683     if (track) track->AddFlow(v1,v2,v3,v4,v5,rp1,rp2,rp3,rp4,rp5,fAfterBurnerPrecision);
 684   }
 685   SetUserModified();
 686 }
 687
 688 //_____________________________________________________________________________
 689 void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5 )
 690 {
 691   //add flow to all tracks wrt the reaction plane angle
 692   for (Int_t i=0; i<fNumberOfTracks; i++)
 693   {
 694     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 695     if (track) track->AddFlow(v1,v2,v3,v4,v5,fMCReactionPlaneAngle, fAfterBurnerPrecision);
 696   }
 697   SetUserModified();
 698 }
 699
 700 //_____________________________________________________________________________
 701 void AliFlowEventSimple::AddV2( TF1* ptDepV2 )
 702 {
 703   //add v2 to all tracks wrt the reaction plane angle
 704   for (Int_t i=0; i<fNumberOfTracks; i++)
 705   {
 706     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 707     if (!track) continue;
 708     Double_t v2 = ptDepV2->Eval(track->Pt());
 709     track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 710   }
 711   SetUserModified();
 712 }
 713
 714 //_____________________________________________________________________________
 715 void AliFlowEventSimple::TagRP( const AliFlowTrackSimpleCuts* cuts )
 716 {
 717   //tag tracks as reference particles (RPs)
 718   for (Int_t i=0; i<fNumberOfTracks; i++)
 719   {
 720     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 721     if (!track) continue;
 722     Bool_t pass=cuts->PassesCuts(track);
 723     Bool_t rpTrack=track->InRPSelection();
 724     if (pass)
 725     {
 726       if (!rpTrack) fNumberOfRPs++; //only increase if not already tagged
 727     }
 728     else
 729     {
 730       if (rpTrack) fNumberOfRPs--; //only decrease if detagging
 731     }
 732     track->SetForRPSelection(pass);
 733   }
 734 }
 735
 736 //_____________________________________________________________________________
 737 void AliFlowEventSimple::TagPOI( const AliFlowTrackSimpleCuts* cuts )
 738 {
 739   //tag tracks as particles of interest (POIs)
 740   for (Int_t i=0; i<fNumberOfTracks; i++)
 741   {
 742     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 743     if (!track) continue;
 744     Bool_t pass=cuts->PassesCuts(track);
 745     track->SetForPOISelection(pass);
 746   }
 747 }
 748
 749 //_____________________________________________________________________________
 750 void AliFlowEventSimple::DefineDeadZone( Double_t etaMin,
 751                                          Double_t etaMax,
 752                                          Double_t phiMin,
 753                                          Double_t phiMax )
 754 {
 755   //mark tracks in given eta-phi region as dead
 756   //by resetting the flow bits
 757   for (Int_t i=0; i<fNumberOfTracks; i++)
 758   {
 759     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 760     Double_t eta = track->Eta();
 761     Double_t phi = track->Phi();
 762     if (eta>etaMin && eta<etaMax && phi>phiMin && phi<phiMax)
 763     {
 764       if (track->InRPSelection()) fNumberOfRPs--;
 765       track->ResetFlowTags();
 766     }
 767   }
 768 }
 769
 770 //_____________________________________________________________________________
 771 Int_t AliFlowEventSimple::CleanUpDeadTracks()
 772 {
 773   //remove tracks that have no flow tags set and cleanup the container
 774   //returns number of cleaned tracks
 775   Int_t ncleaned=0;
 776   for (Int_t i=0; i<fNumberOfTracks; i++)
 777   {
 778     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 779     if (!track) continue;
 780     if (track->IsDead()) {delete track;track=NULL;ncleaned++;}
 781   }
 782   fTrackCollection->Compress(); //clean up empty slots
 783   fNumberOfTracks-=ncleaned; //update number of tracks
 784   return ncleaned;
 785 }
 786
 787 //_____________________________________________________________________________
 788 TF1* AliFlowEventSimple::SimplePtDepV2()
 789 {
 790   //return a standard pt dependent v2 formula, user has to clean up!
 791   return new TF1("StandardPtDepV2","((x<1.0)*(0.05/1.0)*x+(x>=1.0)*0.05)");
 792 }
 793
 794 //_____________________________________________________________________________
 795 TF1* AliFlowEventSimple::SimplePtSpectrum()
 796 {
 797   //return a standard pt spectrum, user has to clean up!
 798   return new TF1("StandardPtSpectrum","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
 799 }
 800
 801 //_____________________________________________________________________________
 802 void AliFlowEventSimple::ClearFast()
 803 {
 804   //clear the counter without deleting allocated objects so they can be reused
 805   fReferenceMultiplicity = 0;
 806   fNumberOfTracks = 0;
 807   fNumberOfRPs = 0;
 808   fMCReactionPlaneAngle = 0.0;
 809   fMCReactionPlaneAngleIsSet = kFALSE;
 810   fAfterBurnerPrecision = 0.001;
 811   fUserModified = kFALSE;
 812 }