[u/mrichter/AliRoot.git] / ANALYSIS / AliFlowAnalysis.cxx

#include "AliFlowAnalysis.h"
//________________________________
///////////////////////////////////////////////////////////
//
// class AliFlowAnalysis
//
// Flow Analysis
//
//
// S.Radomski@gsi.de
// Piotr.Skowronski@cern.ch
//
///////////////////////////////////////////////////////////
/*********************************************************/

#include <TString.h>
#include <TParticle.h>

#include <AliStack.h>
#include <AliAOD.h>
#include <AliVAODParticle.h>
#include <AliAODParticleCut.h>

#include <AliESDtrack.h>
#include <AliESD.h>

ClassImp(AliFlowAnalysis)

AliFlowAnalysis::AliFlowAnalysis():
 fPartCut(0x0)
{
 //ctor
}
/*********************************************************/

AliFlowAnalysis::~AliFlowAnalysis()
{
 //dtor
  delete fPartCut;
}
/*********************************************************/

Int_t AliFlowAnalysis::Init()
{
  //Initilizes anaysis
  Info("Init","");
  return 0; 
}
/*********************************************************/

Int_t AliFlowAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
{
 
  Info("ProcessEvent","Sim AOD address %#x",aodsim);
  Double_t psi = 0, v2 = 0;
  if (aodrec)
   {
     GetFlow(aodrec,v2,psi);
     Info("ProcessEvent","Reconstructed Event: Event plane is %f, V2 is %f",psi,v2);
   }  

  if (aodsim)
   {
     GetFlow(aodsim,v2,psi);
     Info("ProcessEvent","Simulated Event: Event plane is %f, V2 is %f",psi,v2);
   }  
  
  return 0;
  
}
/*********************************************************/

Int_t AliFlowAnalysis::Finish()
{
  //Finish analysis and writes results
  Info("Init","Finish");
  return 0;
}
/*********************************************************/

Double_t AliFlowAnalysis::GetEventPlane(AliAOD* aod)
{
  //returns event plane in degrees
  if (aod == 0x0)
   {
     Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
     return -1.0;
   }

  Double_t psi;
  Int_t mult = aod->GetNumberOfParticles();
  
  Double_t ssin = 0, scos = 0;

  for (Int_t i=0; i<mult; i++) 
   {
     AliVAODParticle* aodtrack = aod->GetParticle(i);
     if (aodtrack == 0x0)
      {
        Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
        continue;
      }
     
     if (fPartCut)
      if (fPartCut->Rejected(aodtrack))
        continue;

     Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px()); 
     
     ssin += TMath::Sin( 2.0 * phi );
     scos += TMath::Cos( 2.0 * phi );
   }

  psi = atan2 (ssin, scos) / 2.0;
  psi = psi * 180. / TMath::Pi(); 
  
  return psi;

}
/*********************************************************/

void AliFlowAnalysis::GetFlow(AliAOD* aod,Double_t& v2,Double_t& psi)
{
//returns flow parameters: v2 and event plane
  if (aod == 0x0)
   {
     Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
     return;
   }
   
  psi = GetEventPlane(aod);
  Int_t mult = aod->GetNumberOfParticles();
  
  Double_t ssin = 0, scos = 0;

  for (Int_t i=0; i<mult; i++) 
   {
     AliVAODParticle* aodtrack = aod->GetParticle(i);
     if (aodtrack == 0x0)
      {
        Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
        continue;
      }
     if (fPartCut)
      if (fPartCut->Rejected(aodtrack))
        continue;
      
     Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px()); 
     ssin += TMath::Sin( 2.0 * (phi - psi));
     scos += TMath::Cos( 2.0 * (phi - psi));
   }
   
  v2 = TMath::Hypot(ssin,scos);
}


/*********************************************************/

Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd)
{
  //returns event plane
  if (esd == 0x0)
   {
     ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
     return -1.0;
   }

  Double_t psi;
  Int_t mult = esd->GetNumberOfTracks();
  
  Double_t ssin = 0, scos = 0;

  for (Int_t i=0; i<mult; i++) 
   {
     AliESDtrack* esdtrack = esd->GetTrack(i);
     if (esdtrack == 0x0)
      {
        ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
        continue;
      }
      
     Double_t mom[3];//momentum
     esdtrack->GetPxPyPz(mom); 
     Double_t phi = TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]); 
     
     ssin += TMath::Sin( 2.0 * phi );
     scos += TMath::Cos( 2.0 * phi );
   }

  psi = atan2 (ssin, scos) / 2.0;
  psi = psi * 180. / TMath::Pi(); 
  
  return psi;

}
/*********************************************************/

void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi)
{
//returns flow parameters: v2 and event plane
  if (esd == 0x0)
   {
     ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
     return;
   }
   
  psi = GetEventPlane(esd);
  Int_t mult = esd->GetNumberOfTracks();
  
  Double_t ssin = 0, scos = 0;

  for (Int_t i=0; i<mult; i++) 
   {
     AliESDtrack* esdtrack = esd->GetTrack(i);
     if (esdtrack == 0x0)
      {
        ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
        continue;
      }
      
     Double_t mom[3];//momentum
     esdtrack->GetPxPyPz(mom); 
     Double_t phi = TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]); 
     
     ssin += TMath::Sin( 2.0 * (phi - psi));
     scos += TMath::Cos( 2.0 * (phi - psi));
   }
   
  v2 = TMath::Hypot(ssin,scos);
}
/*********************************************************/
Commit	Line	Data
b26900d0	1	#include "AliFlowAnalysis.h"
	2	//________________________________
	3	///////////////////////////////////////////////////////////
	4	//
	5	// class AliFlowAnalysis
	6	//
	7	// Flow Analysis
	8	//
	9	//
	10	// S.Radomski@gsi.de
	11	// Piotr.Skowronski@cern.ch
	12	//
	13	///////////////////////////////////////////////////////////
	14	/*********************************************************/
	15
	16	#include <TString.h>
	17	#include <TParticle.h>
	18
	19	#include <AliStack.h>
a5556ea5	20	#include <AliAOD.h>
afa8b37b	21	#include <AliVAODParticle.h>
a5556ea5	22	#include <AliAODParticleCut.h>
a5556ea5	23
b26900d0	24	#include <AliESDtrack.h>
	25	#include <AliESD.h>
	26
	27	ClassImp(AliFlowAnalysis)
	28
a5556ea5	29	AliFlowAnalysis::AliFlowAnalysis():
	30	fPartCut(0x0)
	31	{
	32	//ctor
	33	}
	34	/*********************************************************/
	35
	36	AliFlowAnalysis::~AliFlowAnalysis()
	37	{
	38	//dtor
	39	delete fPartCut;
	40	}
	41	/*********************************************************/
	42
b26900d0	43	Int_t AliFlowAnalysis::Init()
	44	{
	45	//Initilizes anaysis
	46	Info("Init","");
	47	return 0;
	48	}
b26900d0	49	/*********************************************************/
b26900d0	50
a5556ea5	51	Int_t AliFlowAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
b26900d0	52	{
b26900d0	53
a5556ea5	54	Info("ProcessEvent","Sim AOD address %#x",aodsim);
	55	Double_t psi = 0, v2 = 0;
	56	if (aodrec)
	57	{
	58	GetFlow(aodrec,v2,psi);
	59	Info("ProcessEvent","Reconstructed Event: Event plane is %f, V2 is %f",psi,v2);
	60	}
	61
	62	if (aodsim)
	63	{
	64	GetFlow(aodsim,v2,psi);
	65	Info("ProcessEvent","Simulated Event: Event plane is %f, V2 is %f",psi,v2);
	66	}
	67
b26900d0	68	return 0;
	69
	70	}
	71	/*********************************************************/
	72
	73	Int_t AliFlowAnalysis::Finish()
	74	{
	75	//Finish analysis and writes results
	76	Info("Init","Finish");
	77	return 0;
	78	}
	79	/*********************************************************/
	80
a5556ea5	81	Double_t AliFlowAnalysis::GetEventPlane(AliAOD* aod)
	82	{
	83	//returns event plane in degrees
	84	if (aod == 0x0)
	85	{
	86	Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
	87	return -1.0;
	88	}
	89
	90	Double_t psi;
	91	Int_t mult = aod->GetNumberOfParticles();
	92
	93	Double_t ssin = 0, scos = 0;
	94
	95	for (Int_t i=0; i<mult; i++)
	96	{
afa8b37b	97	AliVAODParticle* aodtrack = aod->GetParticle(i);
a5556ea5	98	if (aodtrack == 0x0)
	99	{
	100	Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
	101	continue;
	102	}
	103
	104	if (fPartCut)
cea0a066	105	if (fPartCut->Rejected(aodtrack))
a5556ea5	106	continue;
	107
	108	Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px());
	109
	110	ssin += TMath::Sin( 2.0 * phi );
	111	scos += TMath::Cos( 2.0 * phi );
	112	}
	113
	114	psi = atan2 (ssin, scos) / 2.0;
	115	psi = psi * 180. / TMath::Pi();
	116
	117	return psi;
	118
	119	}
	120	/*********************************************************/
	121
	122	void AliFlowAnalysis::GetFlow(AliAOD* aod,Double_t& v2,Double_t& psi)
	123	{
	124	//returns flow parameters: v2 and event plane
	125	if (aod == 0x0)
	126	{
	127	Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
	128	return;
	129	}
	130
	131	psi = GetEventPlane(aod);
	132	Int_t mult = aod->GetNumberOfParticles();
	133
	134	Double_t ssin = 0, scos = 0;
	135
	136	for (Int_t i=0; i<mult; i++)
	137	{
afa8b37b	138	AliVAODParticle* aodtrack = aod->GetParticle(i);
a5556ea5	139	if (aodtrack == 0x0)
	140	{
	141	Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
	142	continue;
	143	}
	144	if (fPartCut)
cea0a066	145	if (fPartCut->Rejected(aodtrack))
a5556ea5	146	continue;
	147
	148	Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px());
	149	ssin += TMath::Sin( 2.0 * (phi - psi));
	150	scos += TMath::Cos( 2.0 * (phi - psi));
	151	}
	152
	153	v2 = TMath::Hypot(ssin,scos);
	154	}
	155
	156
	157	/*********************************************************/
	158
b26900d0	159	Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd)
	160	{
	161	//returns event plane
e6577ac7	162	if (esd == 0x0)
e6577ac7	163	{
a5556ea5	164	::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
e6577ac7	165	return -1.0;
	166	}
	167
b26900d0	168	Double_t psi;
	169	Int_t mult = esd->GetNumberOfTracks();
	170
	171	Double_t ssin = 0, scos = 0;
	172
	173	for (Int_t i=0; i<mult; i++)
	174	{
	175	AliESDtrack* esdtrack = esd->GetTrack(i);
	176	if (esdtrack == 0x0)
	177	{
a5556ea5	178	::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
b26900d0	179	continue;
	180	}
	181
	182	Double_t mom[3];//momentum
	183	esdtrack->GetPxPyPz(mom);
	184	Double_t phi = TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
	185
	186	ssin += TMath::Sin( 2.0 * phi );
	187	scos += TMath::Cos( 2.0 * phi );
	188	}
	189
	190	psi = atan2 (ssin, scos) / 2.0;
	191	psi = psi * 180. / TMath::Pi();
	192
	193	return psi;
	194
	195	}
a5556ea5	196	/*********************************************************/
e6577ac7	197
	198	void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi)
	199	{
	200	//returns flow parameters: v2 and event plane
	201	if (esd == 0x0)
	202	{
a5556ea5	203	::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
e6577ac7	204	return;
	205	}
	206
	207	psi = GetEventPlane(esd);
	208	Int_t mult = esd->GetNumberOfTracks();
	209
	210	Double_t ssin = 0, scos = 0;
	211
	212	for (Int_t i=0; i<mult; i++)
	213	{
	214	AliESDtrack* esdtrack = esd->GetTrack(i);
	215	if (esdtrack == 0x0)
	216	{
a5556ea5	217	::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
e6577ac7	218	continue;
	219	}
	220
	221	Double_t mom[3];//momentum
	222	esdtrack->GetPxPyPz(mom);
	223	Double_t phi = TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
	224
	225	ssin += TMath::Sin( 2.0 * (phi - psi));
	226	scos += TMath::Cos( 2.0 * (phi - psi));
	227	}
	228
	229	v2 = TMath::Hypot(ssin,scos);
	230	}
a5556ea5	231	/*********************************************************/