[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimple.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include "Riostream.h"
#include "TObjArray.h"
#include "TFile.h"
#include "TList.h"
#include "TMath.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TBrowser.h"
#include "AliFlowVector.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowEventSimple.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),
  fNumberOfTracks(0),
  fEventNSelTracksRP(0),
  fMCReactionPlaneAngle(0.),
  fNumberOfTracksWrap(NULL),
  fEventNSelTracksRPWrap(NULL),
  fMCReactionPlaneAngleWrap(NULL)
{
  cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
}

//-----------------------------------------------------------------------

AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
  fTrackCollection(NULL),
  fNumberOfTracks(0),
  fEventNSelTracksRP(0),
  fMCReactionPlaneAngle(0.),
  fNumberOfTracksWrap(NULL),
  fEventNSelTracksRPWrap(NULL),
  fMCReactionPlaneAngleWrap(NULL)
{
  //constructor 
  fTrackCollection =  new TObjArray(aLenght) ;
}

//-----------------------------------------------------------------------

AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
  TObject(),
  fTrackCollection(anEvent.fTrackCollection),
  fNumberOfTracks(anEvent.fNumberOfTracks),
  fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
  fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
  fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
  fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
  fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
{
  //copy constructor 
}

//-----------------------------------------------------------------------

AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
{
  *fTrackCollection = *anEvent.fTrackCollection ;
  fNumberOfTracks = anEvent.fNumberOfTracks;
  fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
  fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
  fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
  fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
  fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;

  return *this;
}

//----------------------------------------------------------------------- 

AliFlowEventSimple::~AliFlowEventSimple()
{
  //destructor
  if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
  if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
  if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
  if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
}

//----------------------------------------------------------------------- 

AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
{
  //get track i from collection
  AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
  return pTrack;
}

//-----------------------------------------------------------------------   
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)  
  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.;
  
  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 
  
  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)
    {
     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)
    {
     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++)                               
  {
   pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
   if(pTrack)
   {
    if(pTrack->InRPSelection()) 
    {
     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->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);

  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) 
{
  //for eta subevents
  
  // 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.;
  
  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 
  
  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)
    {
     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)
    {
     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++)                               
  {
   pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
   if(pTrack)
   {
    if(pTrack->InRPSelection())
    {
     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->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);

  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, fEventNSelTracksRP, fMCReactionPlaneAngle );

  if (fTrackCollection) {  
    fTrackCollection->Print(option);
  }
  else {
    printf( "Empty track collection \n");
  }
}

//----------------------------------------------------------------------- 
 void AliFlowEventSimple::Browse(TBrowser *b)
{
  if (!b) return;
  if (!fNumberOfTracksWrap) {
    fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
    b->Add(fNumberOfTracksWrap);
  }
  if (!fEventNSelTracksRPWrap) {
    fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
    b->Add(fEventNSelTracksRPWrap);
  }
  if (!fMCReactionPlaneAngleWrap) {
    fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
    b->Add( fMCReactionPlaneAngleWrap);
  }
  if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
}
Commit	Line	Data
f1d945a1	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
f1d945a1	16	#include "Riostream.h"
f1d945a1	17	#include "TObjArray.h"
26c4cbb9	18	#include "TFile.h"
03a02aca	19	#include "TList.h"
f1d945a1	20	#include "TMath.h"
26c4cbb9	21	#include "TH1F.h"
	22	#include "TH1D.h"
	23	#include "TProfile.h"
c076fda8	24	#include "TBrowser.h"
f1d945a1	25	#include "AliFlowVector.h"
	26	#include "AliFlowTrackSimple.h"
	27	#include "AliFlowEventSimple.h"
	28
26c4cbb9	29	/**************************************
	30	* AliFlowEventSimple: A simple event *
	31	* for flow analysis *
	32	* *
	33	* authors: Naomi van der Kolk *
	34	* (kolk@nikhef.nl) *
	35	* Ante Bilandzic *
	36	* (anteb@nikhef.nl) *
	37	* ***********************************/
	38
f1d945a1	39	ClassImp(AliFlowEventSimple)
	40
	41	//-----------------------------------------------------------------------
	42
46bec39c	43	AliFlowEventSimple::AliFlowEventSimple():
	44	fTrackCollection(NULL),
	45	fNumberOfTracks(0),
1918addd	46	fEventNSelTracksRP(0),
7183fe85	47	fMCReactionPlaneAngle(0.),
c076fda8	48	fNumberOfTracksWrap(NULL),
1918addd	49	fEventNSelTracksRPWrap(NULL),
a12990bb	50	fMCReactionPlaneAngleWrap(NULL)
46bec39c	51	{
	52	cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
	53	}
	54
	55	//-----------------------------------------------------------------------
	56
	57	AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
	58	fTrackCollection(NULL),
	59	fNumberOfTracks(0),
1918addd	60	fEventNSelTracksRP(0),
7183fe85	61	fMCReactionPlaneAngle(0.),
c076fda8	62	fNumberOfTracksWrap(NULL),
1918addd	63	fEventNSelTracksRPWrap(NULL),
a12990bb	64	fMCReactionPlaneAngleWrap(NULL)
f1d945a1	65	{
f1d945a1	66	//constructor
e35ddff0	67	fTrackCollection = new TObjArray(aLenght) ;
f1d945a1	68	}
	69
	70	//-----------------------------------------------------------------------
	71
e35ddff0	72	AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
bc6b015e	73	TObject(),
e35ddff0	74	fTrackCollection(anEvent.fTrackCollection),
e35ddff0	75	fNumberOfTracks(anEvent.fNumberOfTracks),
1918addd	76	fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
a12990bb	77	fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
c076fda8	78	fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
1918addd	79	fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
a12990bb	80	fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
f1d945a1	81	{
f1d945a1	82	//copy constructor
f1d945a1	83	}
	84
	85	//-----------------------------------------------------------------------
	86
e35ddff0	87	AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
f1d945a1	88	{
26c4cbb9	89	fTrackCollection = anEvent.fTrackCollection ;
e35ddff0	90	fNumberOfTracks = anEvent.fNumberOfTracks;
1918addd	91	fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
a12990bb	92	fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
c076fda8	93	fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
1918addd	94	fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
a12990bb	95	fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
c076fda8	96
f1d945a1	97	return *this;
f1d945a1	98	}
f1d945a1	99
f1d945a1	100	//-----------------------------------------------------------------------
	101
	102	AliFlowEventSimple::~AliFlowEventSimple()
	103	{
	104	//destructor
c076fda8	105	if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
c076fda8	106	if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
1918addd	107	if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
a12990bb	108	if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
f1d945a1	109	}
	110
	111	//-----------------------------------------------------------------------
	112
	113	AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
	114	{
	115	//get track i from collection
e35ddff0	116	AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
e35ddff0	117	return pTrack;
f1d945a1	118	}
	119
	120	//-----------------------------------------------------------------------
03a02aca	121	AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
f1d945a1	122	{
ae733b3b	123	// calculate Q-vector in harmonic n without weights (default harmonic n=2)
e35ddff0	124	Double_t dQX = 0.;
	125	Double_t dQY = 0.;
	126	AliFlowVector vQ;
	127	vQ.Set(0.,0.);
9825d4a9	128
9825d4a9	129	Int_t iOrder = n;
ae733b3b	130	Double_t iUsedTracks = 0;
26c4cbb9	131	Double_t dPhi=0.;
ae733b3b	132	Double_t dPt=0.;
ae733b3b	133	Double_t dEta=0.;
26c4cbb9	134
	135	AliFlowTrackSimple* pTrack = NULL;
	136
	137	Int_t nBinsPhi=0;
	138	Double_t dBinWidthPt=0.;
ae733b3b	139	Double_t dPtMin=0.;
26c4cbb9	140	Double_t dBinWidthEta=0.;
ae733b3b	141	Double_t dEtaMin=0.;
26c4cbb9	142
ae733b3b	143	Double_t wPhi=1.; // weight Phi
	144	Double_t wPt=1.; // weight Pt
	145	Double_t wEta=1.; // weight Eta
26c4cbb9	146
03a02aca	147	TH1F *phiWeights = NULL;
ae733b3b	148	TH1D *ptWeights = NULL;
03a02aca	149	TH1D *etaWeights = NULL;
03a02aca	150
ae733b3b	151	Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 2)
	152	Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 3)
	153	Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 4)
	154	Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 5)
	155	Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 6)
	156	Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 7)
	157	Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 8)
	158
03a02aca	159	if(weightsList)
26c4cbb9	160	{
03a02aca	161	if(usePhiWeights)
26c4cbb9	162	{
03a02aca	163	phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
	164	if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
	165	}
	166	if(usePtWeights)
26c4cbb9	167	{
03a02aca	168	ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
	169	if(ptWeights)
	170	{
	171	dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
ae733b3b	172	dPtMin = (ptWeights->GetXaxis())->GetXmin();
03a02aca	173	}
	174	}
	175	if(useEtaWeights)
26c4cbb9	176	{
03a02aca	177	etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
	178	if(etaWeights)
	179	{
	180	dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
ae733b3b	181	dEtaMin = (etaWeights->GetXaxis())->GetXmin();
03a02aca	182	}
	183	}
	184	} // end of if(weightsList)
26c4cbb9	185
03a02aca	186	// loop over tracks
26c4cbb9	187	for(Int_t i=0;i<fNumberOfTracks;i++)
	188	{
	189	pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i);
	190	if(pTrack)
	191	{
1918addd	192	if(pTrack->InRPSelection())
f1d945a1	193	{
26c4cbb9	194	dPhi = pTrack->Phi();
	195	dPt = pTrack->Pt();
	196	dEta = pTrack->Eta();
ae733b3b	197
	198	// determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
	199	if(phiWeights && nBinsPhi)
26c4cbb9	200	{
ae733b3b	201	wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
26c4cbb9	202	}
ae733b3b	203	// determine v'(pt) weight:
ae733b3b	204	if(ptWeights && dBinWidthPt)
26c4cbb9	205	{
ae733b3b	206	wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
26c4cbb9	207	}
ae733b3b	208	// determine v'(eta) weight:
ae733b3b	209	if(etaWeights && dBinWidthEta)
26c4cbb9	210	{
ae733b3b	211	wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
03a02aca	212	}
ae733b3b	213
ae733b3b	214	// building up the weighted Q-vector:
26c4cbb9	215	dQX += wPhiwPtwEtaTMath::Cos(iOrderdPhi);
26c4cbb9	216	dQY += wPhiwPtwEtaTMath::Sin(iOrderdPhi);
ae733b3b	217
	218	// weighted multiplicity:
	219	iUsedTracks+=wPhiwPtwEta;
	220
	221	// weights raised to various powers are summed up:
	222	dSumOfWeightsToPower2+=pow(wPhiwPtwEta, 2);
	223	dSumOfWeightsToPower3+=pow(wPhiwPtwEta, 3);
	224	dSumOfWeightsToPower4+=pow(wPhiwPtwEta, 4);
	225	dSumOfWeightsToPower5+=pow(wPhiwPtwEta, 5);
	226	dSumOfWeightsToPower6+=pow(wPhiwPtwEta, 6);
	227	dSumOfWeightsToPower7+=pow(wPhiwPtwEta, 7);
	228	dSumOfWeightsToPower8+=pow(wPhiwPtwEta, 8);
03a02aca	229
1918addd	230	} // end of if (pTrack->InRPSelection())
ae733b3b	231	} // end of if (pTrack)
26c4cbb9	232	else {cerr << "no particle!!!"<<endl;}
ae733b3b	233	} // loop over particles
26c4cbb9	234
e35ddff0	235	vQ.Set(dQX,dQY);
e35ddff0	236	vQ.SetMult(iUsedTracks);
ae733b3b	237	vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
	238	vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
	239	vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
	240	vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
	241	vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
	242	vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
	243	vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
26c4cbb9	244
e35ddff0	245	return vQ;
f1d945a1	246
5fef318d	247	}
5fef318d	248
395fadba	249	//-----------------------------------------------------------------------
	250	AliFlowVector AliFlowEventSimple::GetQsub(Double_t etaMin, Double_t etaMax, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
	251	{
	252	//for eta subevents
	253
	254	// calculate Q-vector in harmonic n without weights (default harmonic n=2)
	255	Double_t dQX = 0.;
	256	Double_t dQY = 0.;
	257	AliFlowVector vQ;
	258	vQ.Set(0.,0.);
	259
	260	Int_t iOrder = n;
	261	Double_t iUsedTracks = 0;
	262	Double_t dPhi = 0.;
	263	Double_t dPt = 0.;
	264	Double_t dEta = 0.;
	265
	266	AliFlowTrackSimple* pTrack = NULL;
	267
	268	Int_t nBinsPhi = 0;
	269	Double_t dBinWidthPt = 0.;
	270	Double_t dPtMin = 0.;
	271	Double_t dBinWidthEta= 0.;
	272	Double_t dEtaMin = 0.;
	273
	274	Double_t wPhi = 1.; // weight Phi
	275	Double_t wPt = 1.; // weight Pt
	276	Double_t wEta = 1.; // weight Eta
	277
	278	TH1F *phiWeights = NULL;
	279	TH1D *ptWeights = NULL;
	280	TH1D *etaWeights = NULL;
	281
	282	Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 2)
	283	Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 3)
	284	Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 4)
	285	Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 5)
	286	Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 6)
	287	Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 7)
	288	Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhiwPtwEta)_i, 8)
	289
	290	if(weightsList)
	291	{
	292	if(usePhiWeights)
	293	{
	294	phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
	295	if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
	296	}
	297	if(usePtWeights)
	298	{
	299	ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
	300	if(ptWeights)
	301	{
	302	dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
	303	dPtMin = (ptWeights->GetXaxis())->GetXmin();
	304	}
	305	}
	306	if(useEtaWeights)
	307	{
	308	etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
	309	if(etaWeights)
	310	{
	311	dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
	312	dEtaMin = (etaWeights->GetXaxis())->GetXmin();
313	}
314	}
315	} // end of if(weightsList)
316
317	// loop over tracks
318	for(Int_t i=0;i<fNumberOfTracks;i++)
319	{
320	pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i);
321	if(pTrack)
322	{
1918addd	323	if(pTrack->InRPSelection())
395fadba	324	{
	325	dPhi = pTrack->Phi();
	326	dPt = pTrack->Pt();
	327	dEta = pTrack->Eta();
	328	if (dEta>etaMin && dEta<etaMax) {
	329	// determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
	330	if(phiWeights && nBinsPhi)
	331	{
	332	wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
	333	}
	334	// determine v'(pt) weight:
	335	if(ptWeights && dBinWidthPt)
	336	{
	337	wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
	338	}
	339	// determine v'(eta) weight:
	340	if(etaWeights && dBinWidthEta)
	341	{
	342	wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
	343	}
	344
	345	// building up the weighted Q-vector:
	346	dQX += wPhiwPtwEtaTMath::Cos(iOrderdPhi);
	347	dQY += wPhiwPtwEtaTMath::Sin(iOrderdPhi);
	348
	349	// weighted multiplicity:
	350	iUsedTracks+=wPhiwPtwEta;
	351
	352	// weights raised to various powers are summed up:
	353	dSumOfWeightsToPower2+=pow(wPhiwPtwEta, 2);
	354	dSumOfWeightsToPower3+=pow(wPhiwPtwEta, 3);
	355	dSumOfWeightsToPower4+=pow(wPhiwPtwEta, 4);
	356	dSumOfWeightsToPower5+=pow(wPhiwPtwEta, 5);
	357	dSumOfWeightsToPower6+=pow(wPhiwPtwEta, 6);
	358	dSumOfWeightsToPower7+=pow(wPhiwPtwEta, 7);
	359	dSumOfWeightsToPower8+=pow(wPhiwPtwEta, 8);
	360	} // end of if dEta in eta range
1918addd	361	} // end of if (pTrack->InRPSelection())
395fadba	362	} // end of if (pTrack)
	363	else {cerr << "no particle!!!"<<endl;}
	364	} // loop over particles
	365
	366	vQ.Set(dQX,dQY);
	367	vQ.SetMult(iUsedTracks);
	368	vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
	369	vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
	370	vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
	371	vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
	372	vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
	373	vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
	374	vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
	375
	376	return vQ;
	377
	378	}
	379
	380
c076fda8	381	//-----------------------------------------------------------------------
	382	void AliFlowEventSimple::Print(Option_t *option) const
	383	{
	384	// -----Print some global quantities for this histogram collection class -------*
	385	// ===============================================
	386	// printf( "TH1.Print Name = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
7183fe85	387	printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f",
1918addd	388	GetName(),fNumberOfTracks, fEventNSelTracksRP, fMCReactionPlaneAngle );
c076fda8	389
	390	if (fTrackCollection) {
	391	fTrackCollection->Print(option);
	392	}
	393	else {
	394	printf( "Empty track collection \n");
	395	}
	396	}
	397
	398	//-----------------------------------------------------------------------
	399	void AliFlowEventSimple::Browse(TBrowser *b)
	400	{
	401	if (!b) return;
	402	if (!fNumberOfTracksWrap) {
	403	fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
	404	b->Add(fNumberOfTracksWrap);
	405	}
1918addd	406	if (!fEventNSelTracksRPWrap) {
	407	fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
	408	b->Add(fEventNSelTracksRPWrap);
c076fda8	409	}
a12990bb	410	if (!fMCReactionPlaneAngleWrap) {
7183fe85	411	fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle", fMCReactionPlaneAngle);
a12990bb	412	b->Add( fMCReactionPlaneAngleWrap);
a12990bb	413	}
c076fda8	414	if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
	415	}
	416
395fadba	417