coverity + some patch for pp
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimple.cxx
index fb34022..285a583 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/*****************************************************************
+  AliFlowEventSimple: A simple event 
+  for flow analysis                  
+                                     
+  origin: Naomi van der Kolk (kolk@nikhef.nl)           
+          Ante Bilandzic     (anteb@nikhef.nl)         
+          Raimond Snellings  (Raimond.Snellings@nikhef.nl)    
+  mods:   Mikolaj Krzewicki  (mikolaj.krzewicki@cern.ch)
+*****************************************************************/
+
 #include "Riostream.h"
 #include "TObjArray.h"
 #include "TFile.h"
 #include "TList.h"
+#include "TTree.h"
+#include "TParticle.h"
 #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 "AliFlowTrackSimpleCuts.h"
 #include "AliFlowEventSimple.h"
+#include "TRandom.h"
 
-/**************************************
- * AliFlowEventSimple: A simple event *
- *  for flow analysis                 * 
- *                                    * 
- * authors: Naomi van der Kolk        *
- *           (kolk@nikhef.nl)         *  
- *          Ante Bilandzic            *
- *           (anteb@nikhef.nl)        * 
- * ***********************************/
 ClassImp(AliFlowEventSimple)
 
 //-----------------------------------------------------------------------
-
 AliFlowEventSimple::AliFlowEventSimple():
   fTrackCollection(NULL),
+  fReferenceMultiplicity(0),
   fNumberOfTracks(0),
-  fEventNSelTracksIntFlow(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
+  fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksIntFlowWrap(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),
-  fEventNSelTracksIntFlow(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
+  fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksIntFlowWrap(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(anEvent.fTrackCollection),
+  TObject(anEvent),
+  fTrackCollection((TObjArray*)(anEvent.fTrackCollection)->Clone()),
+  fReferenceMultiplicity(anEvent.fReferenceMultiplicity),
   fNumberOfTracks(anEvent.fNumberOfTracks),
-  fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow),
+  fNumberOfRPs(anEvent.fNumberOfRPs),
   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
+  fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
+  fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
+  fUserModified(anEvent.fUserModified),
   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
-  fEventNSelTracksIntFlowWrap(anEvent.fEventNSelTracksIntFlowWrap),
+  fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 {
-  //copy constructor 
+  //copy constructor
 }
 
 //-----------------------------------------------------------------------
-
 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
 {
-  *fTrackCollection = *anEvent.fTrackCollection ;
+  //assignment operator
+  if (&anEvent==this) return *this; //check self-assignment
+  if (fTrackCollection) fTrackCollection->Delete();
+  delete fTrackCollection;
+  fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
+  fReferenceMultiplicity = anEvent.fReferenceMultiplicity;
   fNumberOfTracks = anEvent.fNumberOfTracks;
-  fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
+  fNumberOfRPs = anEvent.fNumberOfRPs;
   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
-  fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
-  fEventNSelTracksIntFlowWrap = anEvent.fEventNSelTracksIntFlowWrap;
+  fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
+  fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
+  fUserModified=anEvent.fUserModified;
+  fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
+  fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
-
   return *this;
 }
 
-//----------------------------------------------------------------------- 
-
+//-----------------------------------------------------------------------
 AliFlowEventSimple::~AliFlowEventSimple()
 {
   //destructor
-  if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
-  if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
-  if (fEventNSelTracksIntFlowWrap) delete fEventNSelTracksIntFlowWrap;
-  if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
+  if (fTrackCollection) fTrackCollection->Delete();
+  delete fTrackCollection;
+  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();
+}
 
+//-----------------------------------------------------------------------
 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 {
   //get track i from collection
-  AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
+  if (i>=fNumberOfTracks) return NULL;
+  AliFlowTrackSimple* pTrack = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i)) ;
   return pTrack;
 }
 
-//-----------------------------------------------------------------------   
-AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
+//-----------------------------------------------------------------------
+void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* 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++;
+}
+
+//-----------------------------------------------------------------------
+AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 {
-  // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
+  // calculate Q-vector in harmonic n without weights (default harmonic n=2)
   Double_t dQX = 0.;
   Double_t dQY = 0.;
   AliFlowVector vQ;
   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.;
-  Double_t wPhi=1.; // weight Phi  
-  Double_t wPt=1.;  // weight Pt 
-  Double_t wEta=1.; // weight Eta 
-  
+
+  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
+
   TH1F *phiWeights = NULL;
   TH1D *ptWeights  = NULL;
   TH1D *etaWeights = NULL;
 
-  Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
-  Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
-  Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
-  Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
-  Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
-  Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
-  Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
-
   if(weightsList)
   {
-   if(usePhiWeights)
-   {
-    phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
-    if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
-   }          
-   if(usePtWeights)
-   {
-    ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
-    if(ptWeights)
+    if(usePhiWeights)
+    {
+      phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
+      if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
+    }
+    if(usePtWeights)
     {
-     dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
-     dPtMin = (ptWeights->GetXaxis())->GetXmin();
-    } 
-   }       
-   if(useEtaWeights)
-   {
-    etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
-    if(etaWeights)
+      ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
+      if(ptWeights)
+      {
+        dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
+        dPtMin = (ptWeights->GetXaxis())->GetXmin();
+      }
+    }
+    if(useEtaWeights)
     {
-     dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
-     dEtaMin = (etaWeights->GetXaxis())->GetXmin();
-    } 
-   }          
+      etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
+      if(etaWeights)
+      {
+        dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
+        dEtaMin = (etaWeights->GetXaxis())->GetXmin();
+      }
+    }
   } // end of if(weightsList)
-  
-  // loop over tracks    
-  for(Int_t i=0;i<fNumberOfTracks;i++)                               
+
+  // loop over tracks
+  for(Int_t i=0; i<fNumberOfTracks; i++)
   {
-   pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
-   if(pTrack)
-   {
-    if(pTrack->UseForIntegratedFlow()) 
+    pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
+    if(pTrack)
     {
-     dPhi = pTrack->Phi();
-     dPt  = pTrack->Pt();
-     dEta = pTrack->Eta();
-     
-     // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
-     if(phiWeights && nBinsPhi)
-     {
-      wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
-     }
-     // determine v'(pt) weight:    
-     if(ptWeights && dBinWidthPt)
-     {
-      wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
-     }            
-     // determine v'(eta) weight:    
-     if(etaWeights && dBinWidthEta)
-     {
-      wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
-     } 
-
-     // building up the weighted Q-vector:       
-     dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
-     dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
-    
-     // weighted multiplicity:
-     iUsedTracks+=wPhi*wPt*wEta;
-    
-     // weights raised to various powers are summed up:
-     dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
-     dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
-     dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
-     dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
-     dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
-     dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
-     dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
-     
-    } // end of if (pTrack->UseForIntegratedFlow())
-   } // end of if (pTrack)
-   else {cerr << "no particle!!!"<<endl;}
+      if(pTrack->InRPSelection())
+      {
+        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)
+        {
+          wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
+        }
+        // determine v'(pt) weight:
+        if(ptWeights && dBinWidthPt)
+        {
+          wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
+        }
+        // determine v'(eta) weight:
+        if(etaWeights && dBinWidthEta)
+        {
+          wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
+        }
+
+        // building up the weighted Q-vector:
+        dQX += dWeight*wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
+        dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
+
+        // weighted multiplicity:
+        sumOfWeights += dWeight*wPhi*wPt*wEta;
+
+      } // end of if (pTrack->InRPSelection())
+    } // end of if (pTrack)
+    else
+    {
+      cerr << "no particle!!!"<<endl;
+    }
   } // loop over particles
-    
+
   vQ.Set(dQX,dQY);
-  vQ.SetMult(iUsedTracks);
-  vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
-  vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
-  vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
-  vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
-  vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
-  vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
-  vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
+  vQ.SetMult(sumOfWeights);
 
   return vQ;
-  
+
 }
 
-//-----------------------------------------------------------------------   
-AliFlowVector AliFlowEventSimple::GetQsub(Double_t etaMin, Double_t etaMax, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
+//-----------------------------------------------------------------------
+void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 {
-  //for eta subevents
-  
-  // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
+
+  // calculate Q-vector in harmonic n without weights (default harmonic n=2)
   Double_t dQX = 0.;
   Double_t dQY = 0.;
-  AliFlowVector vQ;
-  vQ.Set(0.,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    nBinsPhi    = 0; 
+
+  Int_t    iNbinsPhiSub0 = 0;
+  Int_t    iNbinsPhiSub1 = 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 
-  
-  TH1F *phiWeights = NULL;
-  TH1D *ptWeights  = NULL;
-  TH1D *etaWeights = NULL;
 
-  Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
-  Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
-  Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
-  Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
-  Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
-  Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
-  Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
+  Double_t dWphi = 1.;  // weight Phi
+  Double_t dWpt  = 1.;  // weight Pt
+  Double_t dWeta = 1.;  // weight Eta
+
+  TH1F* phiWeightsSub0 = NULL;
+  TH1F* phiWeightsSub1 = NULL;
+  TH1D* ptWeights  = NULL;
+  TH1D* etaWeights = NULL;
 
   if(weightsList)
+  {
+    if(usePhiWeights)
     {
-      if(usePhiWeights)
-       {
-    phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
-    if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
-   }          
-   if(usePtWeights)
-   {
-    ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
-    if(ptWeights)
+      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)
     {
-     dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
-     dPtMin = (ptWeights->GetXaxis())->GetXmin();
-    } 
-   }       
-   if(useEtaWeights)
-   {
-    etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
-    if(etaWeights)
+      ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
+      if(ptWeights)
+      {
+        dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
+        dPtMin = (ptWeights->GetXaxis())->GetXmin();
+      }
+    }
+    if(useEtaWeights)
     {
-     dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
-     dEtaMin = (etaWeights->GetXaxis())->GetXmin();
-    } 
-   }          
+      etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
+      if(etaWeights)
+      {
+        dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
+        dEtaMin = (etaWeights->GetXaxis())->GetXmin();
+      }
+    }
   } // end of if(weightsList)
-  
-  // loop over tracks    
-  for(Int_t i=0;i<fNumberOfTracks;i++)                               
+
+  //loop over the two subevents
+  for (Int_t s=0; s<2; s++)
   {
-   pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
-   if(pTrack)
-   {
-    if(pTrack->UseForIntegratedFlow())
+    // loop over tracks
+    for(Int_t i=0; i<fNumberOfTracks; i++)
     {
-     dPhi = pTrack->Phi();
-     dPt  = pTrack->Pt();
-     dEta = pTrack->Eta();
-     if (dEta>etaMin && dEta<etaMax) {
-       // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
-       if(phiWeights && nBinsPhi)
-        {
-          wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
-        }
-       // determine v'(pt) weight:    
-       if(ptWeights && dBinWidthPt)
-        {
-          wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
-        }            
-       // determine v'(eta) weight:    
-       if(etaWeights && dBinWidthEta)
-        {
-          wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
-        } 
-
-       // building up the weighted Q-vector:       
-       dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
-       dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
-    
-       // weighted multiplicity:
-       iUsedTracks+=wPhi*wPt*wEta;
-    
-       // weights raised to various powers are summed up:
-       dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
-       dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
-       dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
-       dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
-       dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
-       dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
-       dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
-     } // end of if dEta in eta range
-    } // end of if (pTrack->UseForIntegratedFlow())
-   } // end of if (pTrack)
-   else {cerr << "no particle!!!"<<endl;}
-  } // loop over particles
-    
-  vQ.Set(dQX,dQY);
-  vQ.SetMult(iUsedTracks);
-  vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
-  vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
-  vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
-  vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
-  vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
-  vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
-  vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
+      pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
+      if(pTrack)
+      {
+        if(pTrack->InRPSelection())
+        {
+          if (pTrack->InSubevent(s))
+          {
+            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)
+           //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)
+            {
+              dWeta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
+            }
+
+            // building up the weighted Q-vector:
+            dQX += dWeight*dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
+            dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
+
+            // weighted multiplicity:
+            sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
+
+          } // end of subevent
+        } // end of if (pTrack->InRPSelection())
+      } // end of if (pTrack)
+      else
+      {
+        cerr << "no particle!!!"<<endl;
+      }
+    } // loop over particles
+    Qarray[s].Set(dQX,dQY);
+    Qarray[s].SetMult(sumOfWeights);
+    //reset
+    sumOfWeights = 0.;
+    dQX = 0.;
+    dQY = 0.;
+  }
 
-  return vQ;
-  
 }
 
 
-//----------------------------------------------------------------------- 
+//-----------------------------------------------------------------------
 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",
-         GetName(),fNumberOfTracks, fEventNSelTracksIntFlow, fMCReactionPlaneAngle );
+  printf( "Class.Print Name = %s, #tracks= %d, Number of RPs= %d, MC EventPlaneAngle= %f\n",
+          GetName(),fNumberOfTracks, fNumberOfRPs, fMCReactionPlaneAngle );
 
-  if (fTrackCollection) {  
+  TString optionstr(option);
+  if (!optionstr.Contains("all")) return;
+  if (fTrackCollection)
+  {
     fTrackCollection->Print(option);
   }
-  else {
+  else
+  {
     printf( "Empty track collection \n");
   }
 }
 
-//----------------------------------------------------------------------- 
- void AliFlowEventSimple::Browse(TBrowser *b)
+//-----------------------------------------------------------------------
+void AliFlowEventSimple::Browse(TBrowser *b)
 {
+  //browse in TBrowser
   if (!b) return;
-  if (!fNumberOfTracksWrap) {
+  if (!fNumberOfTracksWrap)
+  {
     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
     b->Add(fNumberOfTracksWrap);
   }
-  if (!fEventNSelTracksIntFlowWrap) {
-    fEventNSelTracksIntFlowWrap = new TParameter<int>("fEventNSelTracksIntFlow", fEventNSelTracksIntFlow);
-    b->Add(fEventNSelTracksIntFlowWrap);
+  if (!fNumberOfRPsWrap)
+  {
+    fNumberOfRPsWrap = new TParameter<int>("fNumberOfRPs", fNumberOfRPs);
+    b->Add(fNumberOfRPsWrap);
   }
-  if (!fMCReactionPlaneAngleWrap) {
+  if (!fMCReactionPlaneAngleWrap)
+  {
     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
     b->Add( fMCReactionPlaneAngleWrap);
   }
   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
 }
 
-  
+//-----------------------------------------------------------------------
+AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
+                                        const AliFlowTrackSimpleCuts* rpCuts,
+                                        const AliFlowTrackSimpleCuts* poiCuts):
+  fTrackCollection(NULL),
+  fReferenceMultiplicity(0),
+  fNumberOfTracks(0),
+  fNumberOfRPs(0),
+  fMCReactionPlaneAngle(0.),
+  fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
+  fNumberOfTracksWrap(NULL),
+  fNumberOfRPsWrap(NULL),
+  fMCReactionPlaneAngleWrap(NULL)
+{
+  //constructor, fills the event from a TTree of kinematic.root files
+  //applies RP and POI cuts, tags the tracks
+
+  Int_t numberOfInputTracks = inputTree->GetEntries() ;
+  fTrackCollection = new TObjArray(numberOfInputTracks/2);
+
+  TParticle* pParticle = new TParticle();
+  inputTree->SetBranchAddress("Particles",&pParticle);
+
+  Int_t iSelParticlesPOI = 0;
+
+  for (Int_t i=0; i<numberOfInputTracks; i++)
+  {
+    inputTree->GetEntry(i);   //get input particle
+    
+    if (!pParticle) continue; //no particle
+    if (!pParticle->IsPrimary()) continue;
+
+    Bool_t rpOK = rpCuts->PassesCuts(pParticle);
+    Bool_t poiOK = poiCuts->PassesCuts(pParticle);
+    
+    if (rpOK || poiOK)
+    {
+      AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);
+
+      //marking the particles used for int. flow:
+      if(rpOK)
+      {
+        pTrack->SetForRPSelection(kTRUE);
+        fNumberOfRPs++;
+      }
+      //marking the particles used for diff. flow:
+      if(poiOK)
+      {
+        pTrack->SetForPOISelection(kTRUE);
+        iSelParticlesPOI++;
+      }
+      //adding a particles which were used either for int. or diff. flow to the list
+      AddTrack(pTrack);
+    }
+  }//for i
+  delete pParticle;
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::CloneTracks(Int_t n)
+{
+  //clone every track n times to add non-flow
+  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<ntracks; itrack++)
+    {
+      AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
+      if (!track) continue;
+      AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
+    }
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::ResolutionPt(Double_t res)
+{
+  //smear pt of all tracks by gaussian with sigma=res
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    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 )
+{
+  //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);
+  }
+}
+
+//_____________________________________________________________________________
+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(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;
+}