added afterburner in flowevent + various fixes

diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.cxx b/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.cxx
index 688d4086ae3047095c32a89a19bf894034ce3a23..f747a7be87113be1f53b58e9160f32d1238c45b8 100644 (file)
--- a/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.cxx
+++ b/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.cxx
@@ -32,6 +32,7 @@
 #include "TMath.h"
 #include "TH1F.h"
 #include "TH1D.h"
+#include "TF1.h"
 #include "TProfile.h"
 #include "TParameter.h"
 #include "TBrowser.h"
@@ -47,45 +48,68 @@ ClassImp(AliFlowEventSimple)
 AliFlowEventSimple::AliFlowEventSimple():
   fTrackCollection(NULL),
   fNumberOfTracks(0),
-  fEventNSelTracksRP(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
   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 aLength):
+  fTrackCollection(new TObjArray(aLength)),
   fNumberOfTracks(0),
-  fEventNSelTracksRP(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksRPWrap(NULL),
+  fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
   //constructor
-  fTrackCollection =  new TObjArray(aLenght);
+}
+
+//-----------------------------------------------------------------------
+AliFlowEventSimple::AliFlowEventSimple( Int_t nParticles,
+                                        TF1* ptDist,
+                                        Double_t phiMin,
+                                        Double_t phiMax,
+                                        Double_t etaMin,
+                                        Double_t etaMax):
+  fTrackCollection(new TObjArray(nParticles)),
+  fNumberOfTracks(0),
+  fNumberOfRPs(0),
+  fMCReactionPlaneAngle(0.),
+  fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
+  fNumberOfTracksWrap(NULL),
+  fNumberOfRPsWrap(NULL),
+  fMCReactionPlaneAngleWrap(NULL)
+{
+  //ctor. generates nParticles random tracks with given Pt distribution
+  //phi and eta are uniform
+  Generate(nParticles,ptDist,phiMin,phiMax,etaMin,etaMax);
 }
 
 //-----------------------------------------------------------------------
 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
-  TObject(),
-  fTrackCollection(NULL),
+  TObject(anEvent),
+  fTrackCollection((TObjArray*)(anEvent.fTrackCollection)->Clone()),
   fNumberOfTracks(anEvent.fNumberOfTracks),
-  fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
+  fNumberOfRPs(anEvent.fNumberOfRPs),
   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
   fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
+  fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
-  fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
+  fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 {
   //copy constructor
-  fTrackCollection = (TObjArray*)anEvent.Clone(); //deep copy
 }
 
 //-----------------------------------------------------------------------
@@ -93,13 +117,14 @@ AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEv
 {
   //assignment operator
   delete fTrackCollection;
-  fTrackCollection = (TObjArray*)anEvent.Clone(); //deep copy
+  fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
   fNumberOfTracks = anEvent.fNumberOfTracks;
-  fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
+  fNumberOfRPs = anEvent.fNumberOfRPs;
   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
   fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
+  fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
-  fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
+  fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
   return *this;
 }
@@ -111,10 +136,28 @@ AliFlowEventSimple::~AliFlowEventSimple()
   if (fTrackCollection) fTrackCollection->Delete();
   delete fTrackCollection;
   if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
-  if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
+  if (fNumberOfRPsWrap) delete fNumberOfRPsWrap;
   if (fMCReactionPlaneAngleWrap) 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
+  for (Int_t i=0; i<nParticles; i++)
+  {
+    AddTrack(new AliFlowTrackSimple( gRandom->Uniform(phiMin,phiMax),
+                                     gRandom->Uniform(etaMin,etaMax),
+                                     ptDist->GetRandom()));
+  }
+}
+
 //-----------------------------------------------------------------------
 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 {
@@ -128,7 +171,6 @@ AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
 {
   //add a track
-  if (!fTrackCollection) return;
   fTrackCollection->AddLast(track);
   fNumberOfTracks++;
 }
@@ -367,7 +409,7 @@ void AliFlowEventSimple::Print(Option_t *option) const
   //             ===============================================
   //   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 );
+          GetName(),fNumberOfTracks, fNumberOfRPs, fMCReactionPlaneAngle );
 
   if (fTrackCollection)
   {
@@ -388,10 +430,10 @@ void AliFlowEventSimple::Browse(TBrowser *b)
     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)
   {
@@ -407,11 +449,12 @@ AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
                                         const AliFlowTrackSimpleCuts* poiCuts):
   fTrackCollection(NULL),
   fNumberOfTracks(0),
-  fEventNSelTracksRP(0),
+  fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
+  fAfterBurnerPrecision(0.001),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksRPWrap(NULL),
+  fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
   //constructor, fills the event from a TTree of kinematic.root files
@@ -443,7 +486,7 @@ AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
       if(rpOK)
       {
         pTrack->SetForRPSelection(kTRUE);
-        fEventNSelTracksRP++;
+        fNumberOfRPs++;
       }
       //marking the particles used for diff. flow:
       if(poiOK)
@@ -462,13 +505,15 @@ 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++)
+  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()));
     }
   }
 }
@@ -497,6 +542,17 @@ void AliFlowEventSimple::TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, D
   }
 }
 
+//_____________________________________________________________________________
+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);
+  }
+}
+
 //_____________________________________________________________________________
 void AliFlowEventSimple::AddV2(Double_t v2)
 {
@@ -504,7 +560,29 @@ void AliFlowEventSimple::AddV2(Double_t v2)
   for (Int_t i=0; i<fNumberOfTracks; i++)
   {
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
-    if (track) track->AddV2(v2, fMCReactionPlaneAngle, 0.001);
+    if (track) track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+}
+
+//_____________________________________________________________________________
+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);
+  }
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddFlow(Double_t v1, Double_t v2, Double_t v4)
+{
+  //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,v4,fMCReactionPlaneAngle, fAfterBurnerPrecision);
   }
 }
 

diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.h b/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.h
index 81c67fc0a2565c7f10bde6a29cfaa9e4fd04de95..059dc08d0e0c36138b2c5cfbd2bbd6adf11417df 100644 (file)
--- a/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.h
+++ b/PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.h
@@ -17,7 +17,9 @@
 
 #include "TObject.h"
 #include "TParameter.h"
+#include "TMath.h"
 class TTree;
+class TF1;
 class AliFlowVector;
 class AliFlowTrackSimple;
 class AliFlowTrackSimpleCuts;
@@ -27,26 +29,44 @@ class AliFlowEventSimple: public TObject {
  public:
   AliFlowEventSimple();
   AliFlowEventSimple(Int_t aLenght);
+  AliFlowEventSimple( Int_t nParticles,
+                      TF1* ptDist,
+                      Double_t phiMin=0,
+                      Double_t phiMax=TMath::TwoPi(),
+                      Double_t etaMin=-1.0,
+                      Double_t etaMax= 1.0 );
   AliFlowEventSimple(TTree* anInput, const AliFlowTrackSimpleCuts* rpCuts, const AliFlowTrackSimpleCuts* poiCuts);
   AliFlowEventSimple(const AliFlowEventSimple& anEvent);
   AliFlowEventSimple& operator=(const AliFlowEventSimple& anEvent);
   virtual  ~AliFlowEventSimple();
 
+  virtual void Generate( Int_t nParticles,
+                         TF1* ptDist,
+                         Double_t phiMin=0,
+                         Double_t phiMax=TMath::TwoPi(),
+                         Double_t etaMin=-1.0,
+                         Double_t etaMax= 1.0 );
+
   Bool_t  IsFolder() const {return kTRUE;};
   void    Browse(TBrowser *b); 
   void    Print(Option_t* option = "") const;      //method to print stats
   
   Int_t    NumberOfTracks() const                   { return fNumberOfTracks; }
-  Int_t    GetEventNSelTracksRP() const             { return fEventNSelTracksRP; } 
-  void     SetEventNSelTracksRP(Int_t nr)           { fEventNSelTracksRP = nr; }  
+  Int_t    GetEventNSelTracksRP() const             { return fNumberOfRPs; } 
+  void     SetEventNSelTracksRP(Int_t nr)           { fNumberOfRPs = nr; }  
   Double_t GetMCReactionPlaneAngle() const          { return fMCReactionPlaneAngle; }
   void     SetMCReactionPlaneAngle(Double_t fPhiRP) { fMCReactionPlaneAngle=fPhiRP; fMCReactionPlaneAngleIsSet=kTRUE; }
   Bool_t   IsSetMCReactionPlaneAngle() const        { return fMCReactionPlaneAngleIsSet; }
+  void     SetAfterBurnerPrecision(Double_t p)      { fAfterBurnerPrecision=p; }
+  Double_t GetAfterBurnerPrecision() const          { return fAfterBurnerPrecision; }
 
   void ResolutionPt(Double_t res);
   void TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, Double_t etaMinB, Double_t etaMaxB );
   void CloneTracks(Int_t n);
+  void AddV1( Double_t v1 );
   void AddV2( Double_t v2 );
+  void AddV4( Double_t v4 );
+  void AddFlow( Double_t v1, Double_t v2, Double_t v4 );
  
   AliFlowTrackSimple* GetTrack(Int_t i);
   void AddTrack( AliFlowTrackSimple* track ); 
@@ -59,16 +79,17 @@ class AliFlowEventSimple: public TObject {
   void     SetNumberOfTracks(Int_t nr)              { fNumberOfTracks = nr; }
 
  protected:
-  TObjArray*              fTrackCollection;           // collection of tracks
+  TObjArray*              fTrackCollection;           //-> collection of tracks
   Int_t                   fNumberOfTracks;            // number of tracks
-  Int_t                   fEventNSelTracksRP;         // number of tracks that have passed the RP selection
+  Int_t                   fNumberOfRPs;               // number of tracks that have passed the RP selection
   Double_t                fMCReactionPlaneAngle;      // the angle of the reaction plane from the MC truth
   Bool_t                  fMCReactionPlaneAngleIsSet; // did we set it from MC?
+  Double_t                fAfterBurnerPrecision;      // iteration precision in afterburner
   TParameter<Int_t>*      fNumberOfTracksWrap;        //! number of tracks in TBrowser
-  TParameter<Int_t>*      fEventNSelTracksRPWrap;     //! number of tracks that have passed the RP selection in TBrowser
-  TParameter<Double_t>*   fMCReactionPlaneAngleWrap;   //! the angle of the reaction plane from the MC truth in TBrowser
+  TParameter<Int_t>*      fNumberOfRPsWrap;           //! number of tracks that have passed the RP selection in TBrowser
+  TParameter<Double_t>*   fMCReactionPlaneAngleWrap;  //! the angle of the reaction plane from the MC truth in TBrowser
 
-  ClassDef(AliFlowEventSimple,1)                       // simplified event used in flow analysis 
+  ClassDef(AliFlowEventSimple,1)
 };
 
 #endif

diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.cxx b/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.cxx
index 868ce171a43c14a983545e75b76fc7da7d747d69..d7f767a65c46feb89ba421f4447571b63c1fd3f7 100644 (file)
--- a/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.cxx
+++ b/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.cxx
@@ -22,7 +22,7 @@
 // author: N. van der Kolk (kolk@nikhef.nl)
 // mods: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
 
-#include "TNamed.h"
+#include "TObject.h"
 #include "TParticle.h"
 #include "AliFlowTrackSimple.h"
 #include "TRandom.h"
@@ -32,6 +32,7 @@ ClassImp(AliFlowTrackSimple)
 
 //-----------------------------------------------------------------------
 AliFlowTrackSimple::AliFlowTrackSimple():
+  TObject(),
   fEta(0),
   fPt(0),
   fPhi(0),
@@ -43,6 +44,7 @@ AliFlowTrackSimple::AliFlowTrackSimple():
 
 //-----------------------------------------------------------------------
 AliFlowTrackSimple::AliFlowTrackSimple(Double_t phi, Double_t eta, Double_t pt):
+  TObject(),
   fEta(eta),
   fPt(pt),
   fPhi(phi),
@@ -54,6 +56,7 @@ AliFlowTrackSimple::AliFlowTrackSimple(Double_t phi, Double_t eta, Double_t pt):
 
 //-----------------------------------------------------------------------
 AliFlowTrackSimple::AliFlowTrackSimple(const TParticle* p):
+  TObject(),
   fEta(p->Eta()),
   fPt(p->Pt()),
   fPhi(p->Phi()),
@@ -64,9 +67,8 @@ AliFlowTrackSimple::AliFlowTrackSimple(const TParticle* p):
 }
 
 //-----------------------------------------------------------------------
-
 AliFlowTrackSimple::AliFlowTrackSimple(const AliFlowTrackSimple& aTrack):
-  TNamed(),
+  TObject(aTrack),
   fEta(aTrack.fEta),
   fPt(aTrack.fPt),
   fPhi(aTrack.fPhi),
@@ -75,8 +77,15 @@ AliFlowTrackSimple::AliFlowTrackSimple(const AliFlowTrackSimple& aTrack):
 {
   //copy constructor 
 }
+
 //-----------------------------------------------------------------------
+AliFlowTrackSimple* AliFlowTrackSimple::Clone(const char* option) const
+{
+  //clone "constructor"
+  return new AliFlowTrackSimple(*this);
+}
 
+//-----------------------------------------------------------------------
 AliFlowTrackSimple& AliFlowTrackSimple::operator=(const AliFlowTrackSimple& aTrack)
 {
   fEta = aTrack.fEta;
@@ -88,7 +97,6 @@ AliFlowTrackSimple& AliFlowTrackSimple::operator=(const AliFlowTrackSimple& aTra
   return *this;
 }
 
-
 //----------------------------------------------------------------------- 
 AliFlowTrackSimple::~AliFlowTrackSimple()
 {
@@ -104,20 +112,99 @@ void AliFlowTrackSimple::ResolutionPt(Double_t res)
 }
 
 //----------------------------------------------------------------------- 
-void AliFlowTrackSimple::AddV2( Double_t v2, Double_t reactionPlaneAngle, Double_t precisionPhi, Int_t maxNumberOfIterations )
+void AliFlowTrackSimple::AddV1( Double_t v1,
+                                Double_t reactionPlaneAngle,
+                                Double_t precisionPhi,
+                                Int_t maxNumberOfIterations )
+{
+  //afterburner, adds v1, uses Newton-Raphson iteration
+  Double_t phi0=fPhi;
+  Double_t f,fp,phiprev;
+
+  for (Int_t i=0; i<maxNumberOfIterations; i++)
+  {
+    phiprev=fPhi; //store last value for comparison
+    f =  fPhi-phi0+2.0*v1*TMath::Sin(fPhi-reactionPlaneAngle);
+    fp = 1.0+2.0*v1*TMath::Cos(fPhi-reactionPlaneAngle); //first derivative
+    fPhi -= f/fp;
+    if (TMath::AreEqualAbs(phiprev,fPhi,precisionPhi)) break;
+  }
+}
+
+//----------------------------------------------------------------------- 
+void AliFlowTrackSimple::AddV2( Double_t v2,
+                                Double_t reactionPlaneAngle,
+                                Double_t precisionPhi,
+                                Int_t maxNumberOfIterations )
 {
   //afterburner, adds v2, uses Newton-Raphson iteration
   Double_t phi0=fPhi;
-  Double_t f,fp,v2sin,v2cos,phiprev;
+  Double_t f,fp,phiprev;
+
+  for (Int_t i=0; i<maxNumberOfIterations; i++)
+  {
+    phiprev=fPhi; //store last value for comparison
+    f =  fPhi-phi0+v2*TMath::Sin(2.*(fPhi-reactionPlaneAngle));
+    fp = 1.0+2.0*v2*TMath::Cos(2.*(fPhi-reactionPlaneAngle)); //first derivative
+    fPhi -= f/fp;
+    if (TMath::AreEqualAbs(phiprev,fPhi,precisionPhi)) break;
+  }
+}
+
+//----------------------------------------------------------------------- 
+void AliFlowTrackSimple::AddV4( Double_t v4,
+                                Double_t reactionPlaneAngle,
+                                Double_t precisionPhi,
+                                Int_t maxNumberOfIterations )
+{
+  //afterburner, adds v4, uses Newton-Raphson iteration
+  Double_t phi0=fPhi;
+  Double_t f,fp,phiprev;
+
+  for (Int_t i=0; i<maxNumberOfIterations; i++)
+  {
+    phiprev=fPhi; //store last value for comparison
+    f =  fPhi-phi0+0.5*v4*TMath::Sin(4.*(fPhi-reactionPlaneAngle));
+    fp = 1.0+2.0*v4*TMath::Cos(4.*(fPhi-reactionPlaneAngle)); //first derivative
+    fPhi -= f/fp;
+    if (TMath::AreEqualAbs(phiprev,fPhi,precisionPhi)) break;
+  }
+}
+
+//______________________________________________________________________________
+void AliFlowTrackSimple::AddFlow( Double_t v1,
+                                  Double_t v2,
+                                  Double_t v4,
+                                  Double_t reactionPlaneAngle,
+                                  Double_t precisionPhi,
+                                  Int_t maxNumberOfIterations )
+{
+  //afterburner, adds v1,v2,v4 uses Newton-Raphson iteration
+  Double_t phi0=fPhi;
+  Double_t f,fp,phiprev;
 
   for (Int_t i=0; i<maxNumberOfIterations; i++)
   {
     phiprev=fPhi; //store last value for comparison
-    v2sin = v2*TMath::Sin(2.*(fPhi-reactionPlaneAngle));
-    v2cos = v2*TMath::Cos(2.*(fPhi-reactionPlaneAngle));
-    f = fPhi-phi0+v2sin;
-    fp = 1.+2.*v2cos; //first derivative
+    f =  fPhi-phi0
+        +2.0*v1*TMath::Sin(fPhi-reactionPlaneAngle)
+        +    v2*TMath::Sin(2.*(fPhi-reactionPlaneAngle))
+        +0.5*v4*TMath::Sin(4.*(fPhi-reactionPlaneAngle))
+        ;
+    fp =  1.0
+         +2.0*(
+           +v1*TMath::Cos(fPhi-reactionPlaneAngle)
+           +v2*TMath::Cos(2.*(fPhi-reactionPlaneAngle))
+           +v4*TMath::Cos(4.*(fPhi-reactionPlaneAngle))
+         ); //first derivative
     fPhi -= f/fp;
     if (TMath::AreEqualAbs(phiprev,fPhi,precisionPhi)) break;
   }
 }
+
+//______________________________________________________________________________
+void AliFlowTrackSimple::Print( Option_t* option ) const
+{
+  //print stuff
+  printf("Phi: %.3f, Eta: %.3f, Pt: %.3f\n",fPhi,fEta,fPt);
+}

diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.h b/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.h
index e2ff59345763641a6bb77697eead1753091ecc0c..9dc146360862c7454fb3cb348ab78decd5331ce7 100644 (file)
--- a/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.h
+++ b/PWG2/FLOW/AliFlowCommon/AliFlowTrackSimple.h
@@ -5,7 +5,7 @@
 #ifndef ALIFLOWTRACKSIMPLE_H
 #define ALIFLOWTRACKSIMPLE_H
 
-#include "TNamed.h"
+#include "TObject.h"
 #include "TBits.h"
 class TParticle;
 
@@ -14,20 +14,20 @@ class TParticle;
 // author: N. van der Kolk (kolk@nikhef.nl)
 // mods: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
 
-class AliFlowTrackSimple: public TNamed {
+class AliFlowTrackSimple: public TObject {
 
 public:
   AliFlowTrackSimple();
   AliFlowTrackSimple(const TParticle* p);
   AliFlowTrackSimple(const AliFlowTrackSimple& aTrack);
   AliFlowTrackSimple(Double_t phi, Double_t eta, Double_t pt);
-  AliFlowTrackSimple& operator=(const AliFlowTrackSimple& aTrack);
+  virtual AliFlowTrackSimple& operator=(const AliFlowTrackSimple& aTrack);
   virtual  ~AliFlowTrackSimple();
+  virtual AliFlowTrackSimple* Clone(const char* option="") const;
   
   Bool_t  IsFolder() const {return kTRUE;};
   //  void Browse(TBrowser *b); 
-  //  void Print(Option_t* option = "") const;      //method to print stats
- 
+  virtual void Print(Option_t* option = "") const;
 
   Double_t Eta() const; 
   Double_t Pt()  const; 
@@ -46,8 +46,24 @@ public:
   
   void ResolutionPt(Double_t resolution);
 
-  void AddV2( Double_t v2, Double_t reactionPlaneAngle,
-                Double_t precision, Int_t maxNumberOfIterations=100 );
+  void AddV1( Double_t v1,
+              Double_t reactionPlaneAngle,
+              Double_t precision,
+              Int_t maxNumberOfIterations=100 );
+  void AddV2( Double_t v2,
+              Double_t reactionPlaneAngle,
+              Double_t precision,
+              Int_t maxNumberOfIterations=100 );
+  void AddV4( Double_t v4,
+              Double_t reactionPlaneAngle,
+              Double_t precision,
+              Int_t maxNumberOfIterations=100 );
+  void AddFlow( Double_t v1,
+                Double_t v2,
+                Double_t v4,
+                Double_t reactionPlaneAngle,
+                Double_t precision,
+                Int_t maxNumberOfIterations=100 );
     
  private:
   Double_t fEta;    // eta

diff --git a/PWG2/FLOW/AliFlowTasks/AliFlowEvent.cxx b/PWG2/FLOW/AliFlowTasks/AliFlowEvent.cxx
index aa2a05b53accb574065b6c71f66e1272517f8c8a..0cc88cb498f0c76e2e9a89ee45a318b8a23dae36 100644 (file)
--- a/PWG2/FLOW/AliFlowTasks/AliFlowEvent.cxx
+++ b/PWG2/FLOW/AliFlowTasks/AliFlowEvent.cxx
@@ -50,7 +50,6 @@ AliFlowEvent::AliFlowEvent():
 }
 
 //-----------------------------------------------------------------------
-
 AliFlowEvent::AliFlowEvent(const AliFlowEvent& event):
   AliFlowEventSimple(event)
 {
@@ -58,7 +57,6 @@ AliFlowEvent::AliFlowEvent(const AliFlowEvent& event):
 }
 
 //-----------------------------------------------------------------------
-
 AliFlowEvent& AliFlowEvent::operator=(const AliFlowEvent& event)
 {
   //assignment operator
@@ -66,6 +64,15 @@ AliFlowEvent& AliFlowEvent::operator=(const AliFlowEvent& event)
   return *this;
 }
 
+//-----------------------------------------------------------------------
+AliFlowTrack* AliFlowEvent::GetTrack(Int_t i)
+{
+  //get track i from collection
+  if (i>=fNumberOfTracks) return NULL;
+  AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i)) ;
+  return pTrack;
+}
+
 //-----------------------------------------------------------------------
 void AliFlowEvent::SetMCReactionPlaneAngle(const AliMCEvent* mcEvent)
 {
@@ -132,16 +139,13 @@ AliFlowEvent::AliFlowEvent( const AliMCEvent* anInput,
     }
     if (!(rpOK||poiOK)) continue;
 
-    //TODO maybe make a class AliFlowTrack with a constructor from AliVParticle
-    AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
-    pTrack->SetEta(pParticle->Eta());
-    pTrack->SetPhi(pParticle->Phi());
-    pTrack->SetPt(pParticle->Pt());
+    AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
+    pTrack->SetSource(AliFlowTrack::kFromMC);
 
     if (rpOK && rpCFManager)
     {
       pTrack->SetForRPSelection(kTRUE);
-      fEventNSelTracksRP++;
+      fNumberOfRPs++;
     }
     if (poiOK && poiCFManager)
     {
@@ -180,17 +184,15 @@ AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
     }
     if (!(rpOK || poiOK)) continue;
 
-    //make new AliFLowTrackSimple
-    AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
-    pTrack->SetPt(pParticle->Pt() );
-    pTrack->SetEta(pParticle->Eta() );
-    pTrack->SetPhi(pParticle->Phi() );
+    //make new AliFLowTrack
+    AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
+    pTrack->SetSource(AliFlowTrack::kFromESD);
 
     //marking the particles used for int. flow:
     if(rpOK && rpCFManager)
     {
       pTrack->SetForRPSelection(kTRUE);
-      fEventNSelTracksRP++;
+      fNumberOfRPs++;
     }
     //marking the particles used for diff. flow:
     if(poiOK && poiCFManager)
@@ -228,16 +230,14 @@ AliFlowEvent::AliFlowEvent( const AliAODEvent* anInput,
     }
     if (!(rpOK || poiOK)) continue;
 
-    //make new AliFlowTrackSimple
-    AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
-    pTrack->SetPt(pParticle->Pt() );
-    pTrack->SetEta(pParticle->Eta() );
-    pTrack->SetPhi(pParticle->Phi() );
+    //make new AliFlowTrack
+    AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
+    pTrack->SetSource(AliFlowTrack::kFromAOD);
 
     if (rpOK && rpCFManager)
     {
       pTrack->SetForRPSelection(kTRUE);
-      fEventNSelTracksRP++;
+      fNumberOfRPs++;
     }
     if (poiOK && poiCFManager)
     {
@@ -335,8 +335,8 @@ AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 
     if (!(rpOK || poiOK)) continue;
 
-    //make new AliFlowTrackSimple
-    AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
+    //make new AliFlowTrack
+    AliFlowTrack* pTrack = new AliFlowTrack();
     if(anOption == kESDkine)   //take the PID from the MC & the kinematics from the ESD
     {
       pTrack->SetPt(pParticle->Pt() );
@@ -352,7 +352,7 @@ AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 
     if (rpOK && rpCFManager)
     {
-      fEventNSelTracksRP++;
+      fNumberOfRPs++;
       pTrack->SetForRPSelection();
     }
     if (poiOK && poiCFManager) pTrack->SetForPOISelection();

diff --git a/PWG2/FLOW/AliFlowTasks/AliFlowEvent.h b/PWG2/FLOW/AliFlowTasks/AliFlowEvent.h
index 18ec5f0fd1ec9c03b1d6f3c61da75bce501a8ce9..80076b512c10d6069e6e8b8ad767daf5c91ec103 100644 (file)
--- a/PWG2/FLOW/AliFlowTasks/AliFlowEvent.h
+++ b/PWG2/FLOW/AliFlowTasks/AliFlowEvent.h
@@ -13,6 +13,7 @@
 
 class TTree;
 class AliFlowTrackSimpleCuts;
+class AliFlowTrack;
 class AliCFManager;
 class AliMCEvent;
 class AliESDEvent;
@@ -52,6 +53,8 @@ public:
 
   void SetMCReactionPlaneAngle(const AliMCEvent* mcEvent);
 
+  AliFlowTrack* GetTrack( Int_t i );
+
   ClassDef(AliFlowEvent,1)
 };
 

diff --git a/PWG2/FLOW/AliFlowTasks/AliFlowTrack.cxx b/PWG2/FLOW/AliFlowTasks/AliFlowTrack.cxx
index af0a7e5b71aaf7a8e1f2da6af16e6f51044bb12d..7289cfb78d6d00108522f0c3bc10d82476a6f3ff 100644 (file)
--- a/PWG2/FLOW/AliFlowTasks/AliFlowTrack.cxx
+++ b/PWG2/FLOW/AliFlowTasks/AliFlowTrack.cxx
@@ -16,19 +16,26 @@
 /* $Id$ */ 
 
 // AliFlowTrack:
-// A simple track class to the the AliFlowEventSimple for flow analysis
-//
-//
-// author: N. van der Kolk (kolk@nikhef.nl)
-// mods: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
+// A track class for use in AliFlowEvent for flow analysis
+// origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
 
+#include "AliVParticle.h"
 #include "AliFlowTrack.h"
 
 ClassImp(AliFlowTrack)
 
 //-----------------------------------------------------------------------
 AliFlowTrack::AliFlowTrack():
-  AliFlowTrackSimple(0),
+  AliFlowTrackSimple(),
+  fTrackSourceBits(),
+  fFMDmultiplicity(0.)
+{
+  //constructor 
+}
+
+//-----------------------------------------------------------------------
+AliFlowTrack::AliFlowTrack(AliVParticle* p):
+  AliFlowTrackSimple(p->Phi(),p->Eta(),p->Pt()),
   fTrackSourceBits(),
   fFMDmultiplicity(0.)
 {
@@ -44,13 +51,39 @@ AliFlowTrack::AliFlowTrack(const AliFlowTrack& aTrack):
   //copy constructor 
 }
 
+//-----------------------------------------------------------------------
+AliFlowTrack* AliFlowTrack::Clone(const char* option) const
+{
+  //clone "constructor"
+  return new AliFlowTrack(*this);
+}
+
 //-----------------------------------------------------------------------
 AliFlowTrack& AliFlowTrack::operator=(const AliFlowTrack& aTrack)
 {
+  //assignment
   AliFlowTrackSimple::operator=(aTrack);
   fTrackSourceBits = aTrack.fTrackSourceBits;
   fFMDmultiplicity = aTrack.fFMDmultiplicity;
+  return *this;
+}
 
+//-----------------------------------------------------------------------
+AliFlowTrackSimple& AliFlowTrack::operator=(const AliFlowTrackSimple& aTrack)
+{
+  //polymorphic assignment
+  AliFlowTrackSimple::operator=(aTrack);
+  const AliFlowTrack* pft = dynamic_cast<const AliFlowTrack*>(&aTrack);
+  if (pft)
+  {
+    fTrackSourceBits = pft->fTrackSourceBits;
+    fFMDmultiplicity = pft->fFMDmultiplicity;
+  }
+  else
+  {
+    fTrackSourceBits.ResetAllBits();
+    fFMDmultiplicity=0.;
+  }
   return *this;
 }
 

diff --git a/PWG2/FLOW/AliFlowTasks/AliFlowTrack.h b/PWG2/FLOW/AliFlowTasks/AliFlowTrack.h
index 79f309cd38ebb45f6ea82dc08c3001ce1e8fcc9b..a9019586c05236617436dbc83b71f41fae713e1d 100644 (file)
--- a/PWG2/FLOW/AliFlowTasks/AliFlowTrack.h
+++ b/PWG2/FLOW/AliFlowTasks/AliFlowTrack.h
@@ -6,6 +6,7 @@
 #define ALIFLOWTRACK_H
 
 #include "AliFlowTrackSimple.h"
+class AliVParticle;
 
 // AliFlowTrack:
 // A track class to the the AliFlowEvent for flow analysis
@@ -14,14 +15,18 @@
 class AliFlowTrack: public AliFlowTrackSimple {
 
 public:
-  enum trackSource {kFromESD=0,
-                    kFromMC=1,
-                    kFromAOD=2,
-                    kFromTracklet=3};
+  enum trackSource { kFromESD=0,
+                     kFromMC=1,
+                     kFromAOD=2,
+                     kFromTracklet=3,
+                     kFromFMD=4 };
   AliFlowTrack();
+  AliFlowTrack(AliVParticle* p);
   AliFlowTrack& operator=(const AliFlowTrack& aTrack);
+  virtual AliFlowTrackSimple& operator=(const AliFlowTrackSimple& aTrack);
   AliFlowTrack(const AliFlowTrack& aTrack);
   virtual  ~AliFlowTrack();
+  virtual AliFlowTrack* Clone(const char* option="") const;
  
   void SetFMDMultiplicity( const Float_t m ) {fFMDmultiplicity=m;} 
   Float_t GetFMDMultiplicity() const {return fFMDmultiplicity;}
@@ -30,8 +35,6 @@ public:
   Bool_t IsSource( trackSource s ) const
                  { return fTrackSourceBits.TestBitNumber(s); }
 
-
-
 private:
   TBits fTrackSourceBits; //where do i come from?
   Float_t fFMDmultiplicity; //FMD multiplicity