coverity + some patch for pp
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimple.cxx
index 56b8f17..285a583 100644 (file)
 #include "TMath.h"
 #include "TH1F.h"
 #include "TH1D.h"
+#include "TF1.h"
 #include "TProfile.h"
 #include "TParameter.h"
 #include "TBrowser.h"
 #include "AliFlowVector.h"
 #include "AliFlowTrackSimple.h"
-#include "AliFlowEventSimple.h"
 #include "AliFlowTrackSimpleCuts.h"
+#include "AliFlowEventSimple.h"
 #include "TRandom.h"
 
 ClassImp(AliFlowEventSimple)
@@ -46,60 +47,84 @@ ClassImp(AliFlowEventSimple)
 //-----------------------------------------------------------------------
 AliFlowEventSimple::AliFlowEventSimple():
   fTrackCollection(NULL),
+  fReferenceMultiplicity(0),
   fNumberOfTracks(0),
-  fEventNSelTracksRP(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksRPWrap(NULL),
+  fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
 }
 
 //-----------------------------------------------------------------------
-AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
-  fTrackCollection(NULL),
+AliFlowEventSimple::AliFlowEventSimple( Int_t n,
+                                        ConstructionMethod method,
+                                        TF1* ptDist,
+                                        Double_t phiMin,
+                                        Double_t phiMax,
+                                        Double_t etaMin,
+                                        Double_t etaMax):
+  fTrackCollection(new TObjArray(n)),
+  fReferenceMultiplicity(0),
   fNumberOfTracks(0),
-  fEventNSelTracksRP(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksRPWrap(NULL),
+  fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
-  //constructor
-  fTrackCollection =  new TObjArray(aLenght);
+  //ctor
+  // if second argument is set to AliFlowEventSimple::kGenerate
+  // it generates n random tracks with given Pt distribution
+  // (a sane default is provided), phi and eta are uniform
+
+  if (method==kGenerate)
+    Generate(n,ptDist,phiMin,phiMax,etaMin,etaMax);
 }
 
 //-----------------------------------------------------------------------
 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
-  TObject(),
-  fTrackCollection(NULL),
+  TObject(anEvent),
+  fTrackCollection((TObjArray*)(anEvent.fTrackCollection)->Clone()),
+  fReferenceMultiplicity(anEvent.fReferenceMultiplicity),
   fNumberOfTracks(anEvent.fNumberOfTracks),
-  fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
+  fNumberOfRPs(anEvent.fNumberOfRPs),
   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
   fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
+  fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
+  fUserModified(anEvent.fUserModified),
   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
-  fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
+  fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 {
   //copy constructor
-  fTrackCollection = (TObjArray*)anEvent.Clone(); //deep copy
 }
 
 //-----------------------------------------------------------------------
 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
 {
   //assignment operator
+  if (&anEvent==this) return *this; //check self-assignment
+  if (fTrackCollection) fTrackCollection->Delete();
   delete fTrackCollection;
-  fTrackCollection = (TObjArray*)anEvent.Clone(); //deep copy
+  fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
+  fReferenceMultiplicity = anEvent.fReferenceMultiplicity;
   fNumberOfTracks = anEvent.fNumberOfTracks;
-  fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
+  fNumberOfRPs = anEvent.fNumberOfRPs;
   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
   fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
+  fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
+  fUserModified=anEvent.fUserModified;
   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
-  fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
+  fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
   return *this;
 }
@@ -110,9 +135,39 @@ AliFlowEventSimple::~AliFlowEventSimple()
   //destructor
   if (fTrackCollection) fTrackCollection->Delete();
   delete fTrackCollection;
-  if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
-  if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
-  if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
+  delete fNumberOfTracksWrap;
+  delete fNumberOfRPsWrap;
+  delete fMCReactionPlaneAngleWrap;
+}
+
+//-----------------------------------------------------------------------
+void AliFlowEventSimple::Generate(Int_t nParticles,
+                                  TF1* ptDist,
+                                  Double_t phiMin,
+                                  Double_t phiMax,
+                                  Double_t etaMin,
+                                  Double_t etaMax)
+{
+  //generate nParticles random tracks uniform in phi and eta
+  //according to the specified pt distribution
+  if (!ptDist)
+  {
+    static TF1 ptdistribution("ptSpectra","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
+    ptDist=&ptdistribution;
+  }
+
+  for (Int_t i=0; i<nParticles; i++)
+  {
+    AliFlowTrackSimple* track = new AliFlowTrackSimple();
+    track->SetPhi( gRandom->Uniform(phiMin,phiMax) );
+    track->SetEta( gRandom->Uniform(etaMin,etaMax) );
+    track->SetPt( ptDist->GetRandom() );
+    track->SetCharge( (gRandom->Uniform()-0.5<0)?-1:1 );
+    AddTrack(track);
+  }
+  fMCReactionPlaneAngle=gRandom->Uniform(0.0,TMath::TwoPi());
+  fMCReactionPlaneAngleIsSet=kTRUE;
+  SetUserModified();
 }
 
 //-----------------------------------------------------------------------
@@ -127,9 +182,13 @@ AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 //-----------------------------------------------------------------------
 void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
 {
-  //add a track
-  if (!fTrackCollection) return;
-  fTrackCollection->AddLast(track);
+  //add a track, delete the old one if necessary
+  if (fNumberOfTracks < fTrackCollection->GetEntriesFast())
+  {
+    TObject* o = fTrackCollection->At(fNumberOfTracks);
+    if (o) delete o;
+  }
+  fTrackCollection->AddAtAndExpand(track,fNumberOfTracks);
   fNumberOfTracks++;
 }
 
@@ -143,22 +202,23 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
   vQ.Set(0.,0.);
 
   Int_t iOrder = n;
-  Double_t iUsedTracks = 0;
-  Double_t dPhi=0.;
-  Double_t dPt=0.;
-  Double_t dEta=0.;
+  Double_t sumOfWeights = 0.;
+  Double_t dPhi = 0.;
+  Double_t dPt = 0.;
+  Double_t dEta = 0.;
+  Double_t dWeight = 1.;
 
   AliFlowTrackSimple* pTrack = NULL;
 
-  Int_t nBinsPhi=0;
-  Double_t dBinWidthPt=0.;
-  Double_t dPtMin=0.;
-  Double_t dBinWidthEta=0.;
-  Double_t dEtaMin=0.;
+  Int_t nBinsPhi = 0;
+  Double_t dBinWidthPt = 0.;
+  Double_t dPtMin = 0.;
+  Double_t dBinWidthEta = 0.;
+  Double_t dEtaMin = 0.;
 
-  Double_t wPhi=1.; // weight Phi
-  Double_t wPt=1.;  // weight Pt
-  Double_t wEta=1.; // weight Eta
+  Double_t wPhi = 1.; // weight Phi
+  Double_t wPt = 1.;  // weight Pt
+  Double_t wEta = 1.; // weight Eta
 
   TH1F *phiWeights = NULL;
   TH1D *ptWeights  = NULL;
@@ -202,6 +262,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
         dPhi = pTrack->Phi();
         dPt  = pTrack->Pt();
         dEta = pTrack->Eta();
+             dWeight = pTrack->Weight();
 
         // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
         if(phiWeights && nBinsPhi)
@@ -220,11 +281,11 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
         }
 
         // building up the weighted Q-vector:
-        dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
-        dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
+        dQX += dWeight*wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
+        dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
 
         // weighted multiplicity:
-        iUsedTracks+=wPhi*wPt*wEta;
+        sumOfWeights += dWeight*wPhi*wPt*wEta;
 
       } // end of if (pTrack->InRPSelection())
     } // end of if (pTrack)
@@ -235,7 +296,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
   } // loop over particles
 
   vQ.Set(dQX,dQY);
-  vQ.SetMult(iUsedTracks);
+  vQ.SetMult(sumOfWeights);
 
   return vQ;
 
@@ -250,14 +311,16 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
   Double_t dQY = 0.;
 
   Int_t iOrder = n;
-  Double_t iUsedTracks = 0;
+  Double_t sumOfWeights = 0.;
   Double_t dPhi = 0.;
   Double_t dPt  = 0.;
   Double_t dEta = 0.;
+  Double_t dWeight = 1.;
 
   AliFlowTrackSimple* pTrack = NULL;
 
-  Int_t    iNbinsPhi   = 0;
+  Int_t    iNbinsPhiSub0 = 0;
+  Int_t    iNbinsPhiSub1 = 0;
   Double_t dBinWidthPt = 0.;
   Double_t dPtMin      = 0.;
   Double_t dBinWidthEta= 0.;
@@ -267,7 +330,8 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
   Double_t dWpt  = 1.;  // weight Pt
   Double_t dWeta = 1.;  // weight Eta
 
-  TH1F* phiWeights = NULL;
+  TH1F* phiWeightsSub0 = NULL;
+  TH1F* phiWeightsSub1 = NULL;
   TH1D* ptWeights  = NULL;
   TH1D* etaWeights = NULL;
 
@@ -275,10 +339,13 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
   {
     if(usePhiWeights)
     {
-      phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
-      if(phiWeights)
-      {
-        iNbinsPhi = phiWeights->GetNbinsX();
+      phiWeightsSub0 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub0"));
+      if(phiWeightsSub0) {
+       iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
+      }
+      phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
+      if(phiWeightsSub1) {
+       iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
       }
     }
     if(usePtWeights)
@@ -314,20 +381,37 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
         {
           if (pTrack->InSubevent(s))
           {
-            dPhi = pTrack->Phi();
-            dPt  = pTrack->Pt();
-            dEta = pTrack->Eta();
+            dPhi    = pTrack->Phi();
+            dPt     = pTrack->Pt();
+            dEta    = pTrack->Eta();
+           dWeight = pTrack->Weight();
 
             // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
-            if(phiWeights && iNbinsPhi)
-            {
-              dWphi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNbinsPhi/TMath::TwoPi())));
-            }
+           //subevent 0
+           if(s == 0)  { 
+             if(phiWeightsSub0 && iNbinsPhiSub0)  {
+               Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi()));
+               //use the phi value at the center of the bin
+               dPhi  = phiWeightsSub0->GetBinCenter(phiBin);
+               dWphi = phiWeightsSub0->GetBinContent(phiBin);
+             }
+           } 
+           //subevent 1
+           else if (s == 1) { 
+             if(phiWeightsSub1 && iNbinsPhiSub1) {
+               Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi()));
+               //use the phi value at the center of the bin
+               dPhi  = phiWeightsSub1->GetBinCenter(phiBin);
+               dWphi = phiWeightsSub1->GetBinContent(phiBin);
+             } 
+           }
+           
             // determine v'(pt) weight:
             if(ptWeights && dBinWidthPt)
             {
               dWpt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
             }
+
             // determine v'(eta) weight:
             if(etaWeights && dBinWidthEta)
             {
@@ -335,11 +419,11 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
             }
 
             // building up the weighted Q-vector:
-            dQX += dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
-            dQY += dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
+            dQX += dWeight*dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
+            dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
 
             // weighted multiplicity:
-            iUsedTracks+=dWphi*dWpt*dWeta;
+            sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
 
           } // end of subevent
         } // end of if (pTrack->InRPSelection())
@@ -350,9 +434,9 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
       }
     } // loop over particles
     Qarray[s].Set(dQX,dQY);
-    Qarray[s].SetMult(iUsedTracks);
+    Qarray[s].SetMult(sumOfWeights);
     //reset
-    iUsedTracks = 0;
+    sumOfWeights = 0.;
     dQX = 0.;
     dQY = 0.;
   }
@@ -366,9 +450,11 @@ void AliFlowEventSimple::Print(Option_t *option) const
   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
   //             ===============================================
   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
-  printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f\n",
-          GetName(),fNumberOfTracks, fEventNSelTracksRP, fMCReactionPlaneAngle );
+  printf( "Class.Print Name = %s, #tracks= %d, Number of RPs= %d, MC EventPlaneAngle= %f\n",
+          GetName(),fNumberOfTracks, fNumberOfRPs, fMCReactionPlaneAngle );
 
+  TString optionstr(option);
+  if (!optionstr.Contains("all")) return;
   if (fTrackCollection)
   {
     fTrackCollection->Print(option);
@@ -382,16 +468,17 @@ void AliFlowEventSimple::Print(Option_t *option) const
 //-----------------------------------------------------------------------
 void AliFlowEventSimple::Browse(TBrowser *b)
 {
+  //browse in TBrowser
   if (!b) return;
   if (!fNumberOfTracksWrap)
   {
     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
     b->Add(fNumberOfTracksWrap);
   }
-  if (!fEventNSelTracksRPWrap)
+  if (!fNumberOfRPsWrap)
   {
-    fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
-    b->Add(fEventNSelTracksRPWrap);
+    fNumberOfRPsWrap = new TParameter<int>("fNumberOfRPs", fNumberOfRPs);
+    b->Add(fNumberOfRPsWrap);
   }
   if (!fMCReactionPlaneAngleWrap)
   {
@@ -406,12 +493,15 @@ AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
                                         const AliFlowTrackSimpleCuts* rpCuts,
                                         const AliFlowTrackSimpleCuts* poiCuts):
   fTrackCollection(NULL),
+  fReferenceMultiplicity(0),
   fNumberOfTracks(0),
-  fEventNSelTracksRP(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksRPWrap(NULL),
+  fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
   //constructor, fills the event from a TTree of kinematic.root files
@@ -443,7 +533,7 @@ AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
       if(rpOK)
       {
         pTrack->SetForRPSelection(kTRUE);
-        fEventNSelTracksRP++;
+        fNumberOfRPs++;
       }
       //marking the particles used for diff. flow:
       if(poiOK)
@@ -462,15 +552,19 @@ AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
 void AliFlowEventSimple::CloneTracks(Int_t n)
 {
   //clone every track n times to add non-flow
-  for (Int_t i=1; i<n; i++)
+  if (n<=0) return; //no use to clone stuff zero or less times
+  Int_t ntracks = fNumberOfTracks;
+  fTrackCollection->Expand((n+1)*fNumberOfTracks);
+  for (Int_t i=0; i<n; i++)
   {
-    for (Int_t itrack=0; itrack<fNumberOfTracks; itrack++)
+    for (Int_t itrack=0; itrack<ntracks; itrack++)
     {
       AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
       if (!track) continue;
-      AddTrack(new AliFlowTrackSimple(*track));
+      AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
     }
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
@@ -482,15 +576,21 @@ void AliFlowEventSimple::ResolutionPt(Double_t res)
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
     if (track) track->ResolutionPt(res);
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
-void AliFlowEventSimple::TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, Double_t etaMinB, Double_t etaMaxB )
+void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
+                                            Double_t etaMaxA,
+                                            Double_t etaMinB,
+                                            Double_t etaMaxB )
 {
   //Flag two subevents in given eta ranges
   for (Int_t i=0; i<fNumberOfTracks; i++)
   {
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    track->ResetSubEventTags();
     Double_t eta=track->Eta();
     if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
     if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
@@ -498,13 +598,215 @@ void AliFlowEventSimple::TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, D
 }
 
 //_____________________________________________________________________________
-void AliFlowEventSimple::AddFlow(Double_t flow)
+void AliFlowEventSimple::TagSubeventsByCharge()
+{
+  //Flag two subevents in given eta ranges
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    track->ResetSubEventTags();
+    Int_t charge=track->Charge();
+    if (charge<0) track->SetForSubevent(0);
+    if (charge>0) track->SetForSubevent(1);
+  }
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV1( Double_t v1 )
+{
+  //add v2 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV1(v1, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV2( Double_t v2 )
+{
+  //add v2 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV3( Double_t v3 )
+{
+  //add v3 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV3(v3, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV4( Double_t v4 )
+{
+  //add v4 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV4(v4, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV5( Double_t v5 )
+{
+  //add v4 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV5(v5, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5,
+                                  Double_t rp1, Double_t rp2, Double_t rp3, Double_t rp4, Double_t rp5 )
+{
+  //add flow to all tracks wrt the reaction plane angle, for all harmonic separate angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddFlow(v1,v2,v3,v4,v5,rp1,rp2,rp3,rp4,rp5,fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5 )
 {
   //add flow to all tracks wrt the reaction plane angle
   for (Int_t i=0; i<fNumberOfTracks; i++)
   {
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
-    if (track) track->AddFlow(flow, fMCReactionPlaneAngle);
+    if (track) track->AddFlow(v1,v2,v3,v4,v5,fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV2( TF1* ptDepV2 )
+{
+  //add v2 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    Double_t v2 = ptDepV2->Eval(track->Pt());
+    track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::TagRP( const AliFlowTrackSimpleCuts* cuts )
+{
+  //tag tracks as reference particles (RPs)
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    Bool_t pass=cuts->PassesCuts(track);
+    Bool_t rpTrack=track->InRPSelection();
+    if (pass) 
+    {
+      if (!rpTrack) fNumberOfRPs++; //only increase if not already tagged
+    }
+    else
+    {
+      if (rpTrack) fNumberOfRPs--; //only decrease if detagging
+    }
+    track->SetForRPSelection(pass);
+  }
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::TagPOI( const AliFlowTrackSimpleCuts* cuts )
+{
+  //tag tracks as particles of interest (POIs)
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    Bool_t pass=cuts->PassesCuts(track);
+    track->SetForPOISelection(pass);
+  }
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::DefineDeadZone( Double_t etaMin,
+                                         Double_t etaMax,
+                                         Double_t phiMin,
+                                         Double_t phiMax )
+{
+  //mark tracks in given eta-phi region as dead
+  //by resetting the flow bits
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    Double_t eta = track->Eta();
+    Double_t phi = track->Phi();
+    if (eta>etaMin && eta<etaMax && phi>phiMin && phi<phiMax)
+    {
+      if (track->InRPSelection()) fNumberOfRPs--;
+      track->ResetFlowTags();
+    }
   }
 }
 
+//_____________________________________________________________________________
+Int_t AliFlowEventSimple::CleanUpDeadTracks()
+{
+  //remove tracks that have no flow tags set and cleanup the container
+  //returns number of cleaned tracks
+  Int_t ncleaned=0;
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    if (track->IsDead()) {delete track;track=NULL;ncleaned++;}
+  }
+  fTrackCollection->Compress(); //clean up empty slots
+  fNumberOfTracks-=ncleaned; //update number of tracks
+  return ncleaned;
+}
+
+//_____________________________________________________________________________
+TF1* AliFlowEventSimple::SimplePtDepV2()
+{
+  //return a standard pt dependent v2 formula, user has to clean up!
+  return new TF1("StandardPtDepV2","((x<1.0)*(0.05/1.0)*x+(x>=1.0)*0.05)");
+}
+
+//_____________________________________________________________________________
+TF1* AliFlowEventSimple::SimplePtSpectrum()
+{
+  //return a standard pt spectrum, user has to clean up!
+  return new TF1("StandardPtSpectrum","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::ClearFast()
+{
+  //clear the counter without deleting allocated objects so they can be reused
+  fReferenceMultiplicity = 0;
+  fNumberOfTracks = 0;
+  fNumberOfRPs = 0;
+  fMCReactionPlaneAngle = 0.0;
+  fMCReactionPlaneAngleIsSet = kFALSE;
+  fAfterBurnerPrecision = 0.001;
+  fUserModified = kFALSE;
+}