[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / AntiprotonToProton / AliProtonFeedDownAnalysis.cxx

#include <Riostream.h>
#include <TFile.h>
#include <TSystem.h>
#include <TF1.h>
#include <TH2D.h>
#include <TH1D.h>
#include <TH1I.h>
#include <TParticle.h>

#include <AliExternalTrackParam.h>
#include <AliAODEvent.h>
#include <AliESDEvent.h>
//#include <AliLog.h>
#include <AliPID.h>
#include <AliStack.h>
#include <AliCFContainer.h>
#include <AliCFEffGrid.h>
#include <AliCFDataGrid.h>
//#include <AliESDVertex.h>
class AliLog;
class AliESDVertex;

#include "AliProtonFeedDownAnalysis.h"
#include "AliProtonAnalysisBase.h"

ClassImp(AliProtonFeedDownAnalysis)

//____________________________________________________________________//
AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis() : 
  TObject(), fProtonAnalysisBase(0),
  fNBinsY(0), fMinY(0), fMaxY(0),
  fNBinsPt(0), fMinPt(0), fMaxPt(0),
  fProtonContainer(0), fAntiProtonContainer(0), fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
  {
  //Default constructor
 }
//____________________________________________________________________//
/*AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) : 
  TObject(), fProtonAnalysisBase(0),
  fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  fProtonContainer(0), fAntiProtonContainer(0),fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
  {
	//Default constructor
	
	//setting up the containers
	Int_t iBin[2];
	iBin[0] = nbinsY;
	iBin[1] = nbinsPt;
	Double_t *binLimY = new Double_t[nbinsY+1];
	Double_t *binLimPt = new Double_t[nbinsPt+1];
	//values for bin lower bounds
	for(Int_t i = 0; i <= nbinsY; i++) 
		binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
	for(Int_t i = 0; i <= nbinsPt; i++) 
		binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
	
	fProtonContainer = new AliCFContainer("containerProtons","container for protons",4,2,iBin);
	fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
	fProtonContainer->SetBinLimits(1,binLimPt); //pT
	fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",4,2,iBin);
	fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
	fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
	
	
}*/
//____________________________________________________________________//
AliProtonFeedDownAnalysis::~AliProtonFeedDownAnalysis() 
{
	//Default destructor
	if(fProtonAnalysisBase) delete fProtonAnalysisBase;
	if(fProtonContainer) delete fProtonContainer;
	if(fAntiProtonContainer) delete fAntiProtonContainer;
	if(fweightfunction) delete fweightfunction; 
	if(fLambda) delete fLambda;
	if(fLambdaweighted) delete fLambdaweighted;
	if(fAntiLambda) delete fAntiLambda;
	if(fAntiLambdaweighted) delete fAntiLambdaweighted;
}
//____________________________________________________________________//
void AliProtonFeedDownAnalysis::InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) 
{
	//Initializes the histograms
	fNBinsY = nbinsY;
	fMinY = fLowY;
	fMaxY = fHighY;
	fNBinsPt = nbinsPt;
	fMinPt = fLowPt;
	fMaxPt = fHighPt;
	
	
	//setting up the containers
	Int_t iBin[2];
	iBin[0] = nbinsY;
	iBin[1] = nbinsPt;
	Double_t *binLimY = new Double_t[nbinsY+1];
	Double_t *binLimPt = new Double_t[nbinsPt+1];
	//values for bin lower bounds
	for(Int_t i = 0; i <= nbinsY; i++) 
		binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
	for(Int_t i = 0; i <= nbinsPt; i++) 
		binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
	
	fProtonContainer = new AliCFContainer("containerProtons","container for protons",kNSteps,2,iBin);
	fProtonContainer->SetBinLimits(0,binLimY); //rapidity
	fProtonContainer->SetBinLimits(1,binLimPt); //pT
	fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",kNSteps,2,iBin);
	fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
	fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
	/*fLambda=new TH1F("Lambda","Lambda",35,0.5,4.0);
	fLambdaweighted=new TH1F("Lambdaweighted","Lambdaweighted",35,0.5,4.0);
	fAntiLambda=new TH1F("AntiLambda","AntiLambda",35,0.5,4.0);
	fAntiLambdaweighted=new TH1F("AntiLambdaweighted","AntiLambdaweighted",35,0.5,4.0);*/
	fLambda=new TH2F("Lambda","Lambda",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	fLambdaweighted=new TH2F("Lambdaweighted","Lambdaweighted",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	fAntiLambda=new TH2F("AntiLambda","AntiLambda",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	fAntiLambdaweighted=new TH2F("AntiLambdaweighted","AntiLambdaweighted",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	fLambda->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	fLambdaweighted->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	fAntiLambda->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	fAntiLambdaweighted->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	if(fProtonAnalysisBase->GetEtaMode())
	{
		fLambda->GetXaxis()->SetTitle("#eta");
		fLambdaweighted->GetXaxis()->SetTitle("#eta");
		fAntiLambda->GetXaxis()->SetTitle("#eta");
		fAntiLambdaweighted->GetXaxis()->SetTitle("#eta");
	}
	else
	{
		fLambda->GetXaxis()->SetTitle("y");
		fLambdaweighted->GetXaxis()->SetTitle("y");
		fAntiLambda->GetXaxis()->SetTitle("y");
		fAntiLambdaweighted->GetXaxis()->SetTitle("y");
	}	  
}
//_________________________________________________________________________//
void AliProtonFeedDownAnalysis::Analyze(AliESDEvent *esd, const AliESDVertex *vertex,AliStack *stack)
{	
	Int_t nTracks = 0;
	Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;
	Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;
	
	//fHistEvents->Fill(0); //number of analyzed events
	Double_t containerInput[2] ;
	Double_t gPt = 0.0, gP = 0.0;
	nTracks = esd->GetNumberOfTracks();
	Float_t weight;
	Int_t trackflag;
	for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) 
	{
		
		AliESDtrack* track = esd->GetTrack(iTracks);
		AliESDtrack trackTPC;
		Int_t label= track->GetLabel();
	/*	Int_t trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
		if (trackflag!=0)
			weight=GetWeightforProton(label,stack);	
		else
			weight=1.0;	*/
	
		if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) 
		{
			AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
			if(!tpcTrack)
				 continue;
			gPt = tpcTrack->Pt();
			gP = tpcTrack->P();
			if(fProtonAnalysisBase->IsProton(track)) 
			{
				trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
				if (trackflag!=0)
					weight=GetWeightforProton(label,stack);	
				else
					weight=1.0;	
				if(tpcTrack->Charge() > 0) 
				{
					nIdentifiedProtons += 1;
					if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track)) 
					  continue;
					if(!fProtonAnalysisBase->IsAccepted(track)) 
						continue;//track cuts
					if(!fProtonAnalysisBase->IsInPhaseSpace(track)) 
						continue; //track outside the analyzed y-Pt
					nSurvivedProtons += 1;
					if(fProtonAnalysisBase->GetEtaMode()) 
					{
						containerInput[0] = tpcTrack->Eta();
					}
					else 
					{
						//fill the container
						containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
					}
					containerInput[1] = gPt;
					fProtonContainer->Fill(containerInput,kSelected,weight);   
					//Feed-down check
					if (trackflag==1)
					{
						fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
						TParticle *particle  = stack->Particle(label);
						Int_t lmother=particle->GetFirstMother();
						TParticle *mparticle=stack->Particle(lmother);
						Double_t ptlambda= mparticle->Pt();
						Double_t eta_or_y=0.0;
						if(fProtonAnalysisBase->GetEtaMode())
						 	eta_or_y=mparticle->Eta();
						else
							eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
						fLambda->Fill(eta_or_y, ptlambda);
						fLambdaweighted->Fill(eta_or_y, ptlambda,weight);												 
					}	 
					
				}//protons
				else if(tpcTrack->Charge() < 0) 
				{
					nIdentifiedAntiProtons += 1;
					if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track)) 
					  continue;
					if(!fProtonAnalysisBase->IsAccepted(track)) 
						continue;//track cuts
					if(!fProtonAnalysisBase->IsInPhaseSpace(track)) 
						continue; //track outside the analyzed y-Pt
					nSurvivedAntiProtons += 1;
					if(fProtonAnalysisBase->GetEtaMode()) 
					{
						containerInput[0] = tpcTrack->Eta();
					}
					else 
					{
						//fill the container
						containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
					}
					containerInput[1] = gPt;
					fAntiProtonContainer->Fill(containerInput,kSelected,weight);
					//Feed-down check							
					if(trackflag==-1)
					{
						fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
						TParticle *particle  = stack->Particle(label);
						Int_t lmother=particle->GetFirstMother();
						TParticle *mparticle=stack->Particle(lmother);
						Double_t ptlambda= mparticle->Pt();
						Double_t eta_or_y=0.0;
						if(fProtonAnalysisBase->GetEtaMode())
						 	eta_or_y=mparticle->Eta();
						else
							eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
						fAntiLambda->Fill(eta_or_y, ptlambda);
						fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);	
					}						
				}//antiprotons   
			}//proton check
		}//TPC only tracks
		else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal) 
		{
			gPt = track->Pt();
			gP = track->P();
			if(fProtonAnalysisBase->IsProton(track)) 
			{
				trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
				if (trackflag!=0)
					weight=GetWeightforProton(label,stack);	
				else
					weight=1.0;	
				if(track->Charge() > 0) 
				{
					nIdentifiedProtons += 1;
					if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track)) 
					  continue;
					if(!fProtonAnalysisBase->IsAccepted(track)) 
						continue;//track cuts
					if(!fProtonAnalysisBase->IsInPhaseSpace(track)) 
						continue; //track outside the analyzed y-Pt
					nSurvivedProtons += 1;
					if(fProtonAnalysisBase->GetEtaMode()) 
					{
						containerInput[0] = track->Eta();
					}
					else 
					{
						//fill the container
						containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
					}
					containerInput[1] = gPt;
					fProtonContainer->Fill(containerInput,kSelected,weight);  
					 //Feed-down check
					if (trackflag==1)
					{
						fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
						TParticle *particle  = stack->Particle(label);
						Int_t lmother=particle->GetFirstMother();
						TParticle *mparticle=stack->Particle(lmother);
						Double_t ptlambda= mparticle->Pt();
						Double_t eta_or_y=0.0;
						if(fProtonAnalysisBase->GetEtaMode())
						 	eta_or_y=mparticle->Eta();
						else
							eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
						fLambda->Fill(eta_or_y, ptlambda);
						fLambdaweighted->Fill(eta_or_y, ptlambda,weight);		
					}	 
				}//protons
				else if(track->Charge() < 0) 
				{
					nIdentifiedAntiProtons += 1;
					if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track)) 
					  continue;
					if(!fProtonAnalysisBase->IsAccepted(track)) 
						continue;//track cuts
					if(!fProtonAnalysisBase->IsInPhaseSpace(track))
						continue; //track outside the analyzed y-Pt
					nSurvivedAntiProtons += 1;
					if(fProtonAnalysisBase->GetEtaMode()) 
					{
						containerInput[0] = track->Eta();
					}
					else 
					{
						//fill the container
						containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
					}
					containerInput[1] = gPt;
					fAntiProtonContainer->Fill(containerInput,kSelected,weight);   
					if(trackflag==-1)
					{
						fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
						TParticle *particle  = stack->Particle(label);
						Int_t lmother=particle->GetFirstMother();
						TParticle *mparticle=stack->Particle(lmother);
						Double_t ptlambda= mparticle->Pt();
						Double_t eta_or_y=0.0;
						if(fProtonAnalysisBase->GetEtaMode())
						 	eta_or_y=mparticle->Eta();
						else
							eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
						fAntiLambda->Fill(eta_or_y, ptlambda);
						fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
					}	
				}//antiprotons
			}//proton check 
		}//combined tracking
	}//track loop 
	
	if(fProtonAnalysisBase->GetDebugMode())
		Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);

}
//_________________________________________________________________________//
void AliProtonFeedDownAnalysis::Analyze(AliAODEvent *fAOD)
{
  // Analysis from AOD: to be implemented...
  // in the near future.
  fAOD->Print();

}
//_________________________________________________________________________//
void AliProtonFeedDownAnalysis::Analyze(AliStack *stack)
{
	 Double_t containerInput[2] ;
	 Float_t weight;
 	for (Int_t ipart=0; ipart<stack->GetNtrack(); ipart++) 
	{ 
		TParticle *particle  = stack->Particle(ipart);
		Int_t code=particle->GetPdgCode();
		/*if (code==3122)
		{
			fLambda->Fill(particle->Pt());
			fLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));			
		}
		if (code==-3122)
		{
			fAntiLambda->Fill(particle->Pt());
			fAntiLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));			
		}*/
		if (TMath::Abs(code)==2212)
		{
			Int_t trackflag=LambdaIsMother(ipart,stack);//1 mother lambda -1 mother anti lambda 0 mother something else 
			if (trackflag!=0)
				weight=GetWeightforProton(ipart,stack);
			else
				weight=1.0;	
			if(particle->GetUniqueID() == 13) //recject hadronic inter.
				continue; 
			if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) 
				continue;
			if(fProtonAnalysisBase->GetEtaMode()) 
			{
				if((particle->Eta()> fMaxY)||(particle->Eta()< fMinY))
					continue; 
			}	
			else
			{	
				if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz())< fMinY)) 
					continue;
			}		
			if(fProtonAnalysisBase->GetEtaMode()) 
			{
				containerInput[0] =particle->Eta();
			}
			else 
			{
				containerInput[0] = fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz());
			}
			//containerInput[0] = particle->Eta() ;
			containerInput[1] = (Float_t)particle->Pt();
			if (particle->GetPdgCode()==2212)
			{
				fProtonContainer->Fill(containerInput,kAll,weight);
				if(ipart<stack->GetNprimary())
				{
					fProtonContainer->Fill(containerInput,kPrimary,weight);
				}
				else
				{
					if(trackflag==1)
						fProtonContainer->Fill(containerInput,kFromLambda,weight);
				}
			}
			if (particle->GetPdgCode()==-2212)
			{
				fAntiProtonContainer->Fill(containerInput,kAll,weight);
				if(ipart<stack->GetNprimary())
				{
					fAntiProtonContainer->Fill(containerInput,kPrimary,weight);
				}
				else
				{								
					if(trackflag==-1)
						fAntiProtonContainer->Fill(containerInput,kFromLambda,weight);
				}			
			}
		}
	}   
}
//______________________________________________________________________________________________
Int_t AliProtonFeedDownAnalysis::LambdaIsMother(Int_t number, AliStack *stack)
{
	 if(number<0)
	 	return 0;
	TParticle *particle  = stack->Particle(number);
	Int_t motherPDG=0;
	Int_t lmother=-1;
	lmother=particle->GetFirstMother();
	if (lmother<0)		
		return 0;
	TParticle *mparticle=stack->Particle(lmother);
	motherPDG=mparticle->GetPdgCode();										
	if(motherPDG==3122)
		return 1;
	if(motherPDG==-3122)	
		return -1;
	return 0;	
} 
//___________________________________________________________________________________________
Float_t AliProtonFeedDownAnalysis::GetWeightforProton(Int_t number,AliStack *stack)
{
	 if(number<0)
	 	return 1.0;
	TParticle *particle  = stack->Particle(number);
	Int_t lmother=particle->GetFirstMother();
	TParticle *mparticle=stack->Particle(lmother);
	return 	fweightfunction->Eval(mparticle->Pt());
}
//______________________________________________________________________________________________
Float_t AliProtonFeedDownAnalysis::GetWeightforLambda(Float_t pt)
{
	return 	fweightfunction->Eval(pt);
}
Commit	Line	Data
198dfda4	1	#include <Riostream.h>
	2	#include <TFile.h>
	3	#include <TSystem.h>
	4	#include <TF1.h>
	5	#include <TH2D.h>
	6	#include <TH1D.h>
	7	#include <TH1I.h>
	8	#include <TParticle.h>
	9
	10	#include <AliExternalTrackParam.h>
	11	#include <AliAODEvent.h>
	12	#include <AliESDEvent.h>
	13	//#include <AliLog.h>
	14	#include <AliPID.h>
	15	#include <AliStack.h>
	16	#include <AliCFContainer.h>
	17	#include <AliCFEffGrid.h>
	18	#include <AliCFDataGrid.h>
	19	//#include <AliESDVertex.h>
	20	class AliLog;
	21	class AliESDVertex;
	22
	23	#include "AliProtonFeedDownAnalysis.h"
	24	#include "AliProtonAnalysisBase.h"
	25
	26	ClassImp(AliProtonFeedDownAnalysis)
	27
	28	//____________________________________________________________________//
	29	AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis() :
	30	TObject(), fProtonAnalysisBase(0),
	31	fNBinsY(0), fMinY(0), fMaxY(0),
	32	fNBinsPt(0), fMinPt(0), fMaxPt(0),
	33	fProtonContainer(0), fAntiProtonContainer(0), fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
	34	{
	35	//Default constructor
	36	}
	37	//____________________________________________________________________//
c9fcf8c9	38	/*AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
198dfda4	39	TObject(), fProtonAnalysisBase(0),
	40	fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
	41	fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
	42	fProtonContainer(0), fAntiProtonContainer(0),fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)
	43	{
	44	//Default constructor
	45
	46	//setting up the containers
	47	Int_t iBin[2];
	48	iBin[0] = nbinsY;
	49	iBin[1] = nbinsPt;
	50	Double_t *binLimY = new Double_t[nbinsY+1];
	51	Double_t *binLimPt = new Double_t[nbinsPt+1];
	52	//values for bin lower bounds
	53	for(Int_t i = 0; i <= nbinsY; i++)
	54	binLimY[i]=(Double_t)fLowY + (fHighY - fLowY) /nbinsY*(Double_t)i;
	55	for(Int_t i = 0; i <= nbinsPt; i++)
	56	binLimPt[i]=(Double_t)fLowPt + (fHighPt - fLowPt) /nbinsPt*(Double_t)i;
	57
	58	fProtonContainer = new AliCFContainer("containerProtons","container for protons",4,2,iBin);
	59	fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
	60	fProtonContainer->SetBinLimits(1,binLimPt); //pT
	61	fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",4,2,iBin);
	62	fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
	63	fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
	64
c9fcf8c9	65
c9fcf8c9	66	}*/
198dfda4	67	//____________________________________________________________________//
	68	AliProtonFeedDownAnalysis::~AliProtonFeedDownAnalysis()
	69	{
	70	//Default destructor
	71	if(fProtonAnalysisBase) delete fProtonAnalysisBase;
	72	if(fProtonContainer) delete fProtonContainer;
	73	if(fAntiProtonContainer) delete fAntiProtonContainer;
	74	if(fweightfunction) delete fweightfunction;
	75	if(fLambda) delete fLambda;
	76	if(fLambdaweighted) delete fLambdaweighted;
	77	if(fAntiLambda) delete fAntiLambda;
	78	if(fAntiLambdaweighted) delete fAntiLambdaweighted;
	79	}
	80	//____________________________________________________________________//
	81	void AliProtonFeedDownAnalysis::InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt)
	82	{
	83	//Initializes the histograms
	84	fNBinsY = nbinsY;
	85	fMinY = fLowY;
	86	fMaxY = fHighY;
	87	fNBinsPt = nbinsPt;
	88	fMinPt = fLowPt;
	89	fMaxPt = fHighPt;
	90
	91
	92	//setting up the containers
	93	Int_t iBin[2];
	94	iBin[0] = nbinsY;
	95	iBin[1] = nbinsPt;
	96	Double_t *binLimY = new Double_t[nbinsY+1];
	97	Double_t *binLimPt = new Double_t[nbinsPt+1];
	98	//values for bin lower bounds
	99	for(Int_t i = 0; i <= nbinsY; i++)
	100	binLimY[i]=(Double_t)fLowY + (fHighY - fLowY) /nbinsY*(Double_t)i;
	101	for(Int_t i = 0; i <= nbinsPt; i++)
	102	binLimPt[i]=(Double_t)fLowPt + (fHighPt - fLowPt) /nbinsPt*(Double_t)i;
	103
	104	fProtonContainer = new AliCFContainer("containerProtons","container for protons",kNSteps,2,iBin);
	105	fProtonContainer->SetBinLimits(0,binLimY); //rapidity
	106	fProtonContainer->SetBinLimits(1,binLimPt); //pT
	107	fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",kNSteps,2,iBin);
	108	fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
	109	fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
c9fcf8c9	110	/*fLambda=new TH1F("Lambda","Lambda",35,0.5,4.0);
198dfda4	111	fLambdaweighted=new TH1F("Lambdaweighted","Lambdaweighted",35,0.5,4.0);
198dfda4	112	fAntiLambda=new TH1F("AntiLambda","AntiLambda",35,0.5,4.0);
c9fcf8c9	113	fAntiLambdaweighted=new TH1F("AntiLambdaweighted","AntiLambdaweighted",35,0.5,4.0);*/
	114	fLambda=new TH2F("Lambda","Lambda",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	115	fLambdaweighted=new TH2F("Lambdaweighted","Lambdaweighted",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	116	fAntiLambda=new TH2F("AntiLambda","AntiLambda",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	117	fAntiLambdaweighted=new TH2F("AntiLambdaweighted","AntiLambdaweighted",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
	118	fLambda->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	119	fLambdaweighted->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	120	fAntiLambda->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	121	fAntiLambdaweighted->GetYaxis()->SetTitle("P_{T} [GeV/c]");
	122	if(fProtonAnalysisBase->GetEtaMode())
	123	{
	124	fLambda->GetXaxis()->SetTitle("#eta");
	125	fLambdaweighted->GetXaxis()->SetTitle("#eta");
	126	fAntiLambda->GetXaxis()->SetTitle("#eta");
	127	fAntiLambdaweighted->GetXaxis()->SetTitle("#eta");
	128	}
	129	else
	130	{
	131	fLambda->GetXaxis()->SetTitle("y");
	132	fLambdaweighted->GetXaxis()->SetTitle("y");
	133	fAntiLambda->GetXaxis()->SetTitle("y");
	134	fAntiLambdaweighted->GetXaxis()->SetTitle("y");
	135	}
198dfda4	136	}
	137	//_________________________________________________________________________//
	138	void AliProtonFeedDownAnalysis::Analyze(AliESDEvent esd, const AliESDVertex vertex,AliStack *stack)
	139	{
	140	Int_t nTracks = 0;
	141	Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;
	142	Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;
	143
	144	//fHistEvents->Fill(0); //number of analyzed events
	145	Double_t containerInput[2] ;
	146	Double_t gPt = 0.0, gP = 0.0;
	147	nTracks = esd->GetNumberOfTracks();
	148	Float_t weight;
	149	Int_t trackflag;
	150	for(Int_t iTracks = 0; iTracks < nTracks; iTracks++)
	151	{
	152
	153	AliESDtrack* track = esd->GetTrack(iTracks);
	154	AliESDtrack trackTPC;
	155	Int_t label= track->GetLabel();
	156	/* Int_t trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
	157	if (trackflag!=0)
	158	weight=GetWeightforProton(label,stack);
	159	else
	160	weight=1.0; */
	161
	162	if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)\|\|(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid))
	163	{
	164	AliExternalTrackParam tpcTrack = (AliExternalTrackParam )track->GetTPCInnerParam();
	165	if(!tpcTrack)
	166	continue;
	167	gPt = tpcTrack->Pt();
	168	gP = tpcTrack->P();
	169	if(fProtonAnalysisBase->IsProton(track))
	170	{
	171	trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
	172	if (trackflag!=0)
	173	weight=GetWeightforProton(label,stack);
	174	else
	175	weight=1.0;
	176	if(tpcTrack->Charge() > 0)
	177	{
	178	nIdentifiedProtons += 1;
30f87a5b	179	if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
	180	continue;
	181	if(!fProtonAnalysisBase->IsAccepted(track))
198dfda4	182	continue;//track cuts
	183	if(!fProtonAnalysisBase->IsInPhaseSpace(track))
	184	continue; //track outside the analyzed y-Pt
	185	nSurvivedProtons += 1;
	186	if(fProtonAnalysisBase->GetEtaMode())
	187	{
	188	containerInput[0] = tpcTrack->Eta();
	189	}
	190	else
	191	{
	192	//fill the container
	193	containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
	194	}
	195	containerInput[1] = gPt;
	196	fProtonContainer->Fill(containerInput,kSelected,weight);
	197	//Feed-down check
	198	if (trackflag==1)
c9fcf8c9	199	{
	200	fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
	201	TParticle *particle = stack->Particle(label);
	202	Int_t lmother=particle->GetFirstMother();
	203	TParticle *mparticle=stack->Particle(lmother);
	204	Double_t ptlambda= mparticle->Pt();
	205	Double_t eta_or_y=0.0;
	206	if(fProtonAnalysisBase->GetEtaMode())
	207	eta_or_y=mparticle->Eta();
	208	else
	209	eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
	210	fLambda->Fill(eta_or_y, ptlambda);
	211	fLambdaweighted->Fill(eta_or_y, ptlambda,weight);
	212	}
198dfda4	213
	214	}//protons
	215	else if(tpcTrack->Charge() < 0)
	216	{
	217	nIdentifiedAntiProtons += 1;
30f87a5b	218	if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
	219	continue;
	220	if(!fProtonAnalysisBase->IsAccepted(track))
198dfda4	221	continue;//track cuts
	222	if(!fProtonAnalysisBase->IsInPhaseSpace(track))
	223	continue; //track outside the analyzed y-Pt
	224	nSurvivedAntiProtons += 1;
	225	if(fProtonAnalysisBase->GetEtaMode())
	226	{
	227	containerInput[0] = tpcTrack->Eta();
	228	}
	229	else
	230	{
	231	//fill the container
	232	containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());
	233	}
	234	containerInput[1] = gPt;
	235	fAntiProtonContainer->Fill(containerInput,kSelected,weight);
	236	//Feed-down check
	237	if(trackflag==-1)
c9fcf8c9	238	{
	239	fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
	240	TParticle *particle = stack->Particle(label);
	241	Int_t lmother=particle->GetFirstMother();
	242	TParticle *mparticle=stack->Particle(lmother);
	243	Double_t ptlambda= mparticle->Pt();
	244	Double_t eta_or_y=0.0;
	245	if(fProtonAnalysisBase->GetEtaMode())
	246	eta_or_y=mparticle->Eta();
	247	else
	248	eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
	249	fAntiLambda->Fill(eta_or_y, ptlambda);
	250	fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
	251	}
198dfda4	252	}//antiprotons
	253	}//proton check
	254	}//TPC only tracks
	255	else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal)
	256	{
	257	gPt = track->Pt();
	258	gP = track->P();
	259	if(fProtonAnalysisBase->IsProton(track))
	260	{
	261	trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
	262	if (trackflag!=0)
	263	weight=GetWeightforProton(label,stack);
	264	else
	265	weight=1.0;
	266	if(track->Charge() > 0)
	267	{
	268	nIdentifiedProtons += 1;
30f87a5b	269	if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
	270	continue;
	271	if(!fProtonAnalysisBase->IsAccepted(track))
198dfda4	272	continue;//track cuts
	273	if(!fProtonAnalysisBase->IsInPhaseSpace(track))
	274	continue; //track outside the analyzed y-Pt
	275	nSurvivedProtons += 1;
	276	if(fProtonAnalysisBase->GetEtaMode())
	277	{
	278	containerInput[0] = track->Eta();
	279	}
	280	else
	281	{
	282	//fill the container
	283	containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
	284	}
	285	containerInput[1] = gPt;
	286	fProtonContainer->Fill(containerInput,kSelected,weight);
	287	//Feed-down check
	288	if (trackflag==1)
c9fcf8c9	289	{
	290	fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
	291	TParticle *particle = stack->Particle(label);
	292	Int_t lmother=particle->GetFirstMother();
	293	TParticle *mparticle=stack->Particle(lmother);
	294	Double_t ptlambda= mparticle->Pt();
	295	Double_t eta_or_y=0.0;
	296	if(fProtonAnalysisBase->GetEtaMode())
	297	eta_or_y=mparticle->Eta();
	298	else
	299	eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
	300	fLambda->Fill(eta_or_y, ptlambda);
	301	fLambdaweighted->Fill(eta_or_y, ptlambda,weight);
	302	}
198dfda4	303	}//protons
	304	else if(track->Charge() < 0)
	305	{
	306	nIdentifiedAntiProtons += 1;
30f87a5b	307	if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track))
	308	continue;
	309	if(!fProtonAnalysisBase->IsAccepted(track))
198dfda4	310	continue;//track cuts
	311	if(!fProtonAnalysisBase->IsInPhaseSpace(track))
	312	continue; //track outside the analyzed y-Pt
	313	nSurvivedAntiProtons += 1;
	314	if(fProtonAnalysisBase->GetEtaMode())
	315	{
	316	containerInput[0] = track->Eta();
	317	}
	318	else
	319	{
	320	//fill the container
	321	containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());
	322	}
	323	containerInput[1] = gPt;
	324	fAntiProtonContainer->Fill(containerInput,kSelected,weight);
	325	if(trackflag==-1)
c9fcf8c9	326	{
198dfda4	327	fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);
c9fcf8c9	328	TParticle *particle = stack->Particle(label);
	329	Int_t lmother=particle->GetFirstMother();
	330	TParticle *mparticle=stack->Particle(lmother);
	331	Double_t ptlambda= mparticle->Pt();
	332	Double_t eta_or_y=0.0;
	333	if(fProtonAnalysisBase->GetEtaMode())
	334	eta_or_y=mparticle->Eta();
	335	else
	336	eta_or_y=0.5*TMath::Log((mparticle->Energy()+mparticle->Pz())/(mparticle->Energy()-mparticle->Pz()));
	337	fAntiLambda->Fill(eta_or_y, ptlambda);
	338	fAntiLambdaweighted->Fill(eta_or_y, ptlambda,weight);
	339	}
198dfda4	340	}//antiprotons
	341	}//proton check
	342	}//combined tracking
	343	}//track loop
	344
	345	if(fProtonAnalysisBase->GetDebugMode())
	346	Printf("Initial number of tracks: %d \| Identified (anti)protons: %d - %d \| Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
	347
	348	}
	349	//_________________________________________________________________________//
	350	void AliProtonFeedDownAnalysis::Analyze(AliAODEvent *fAOD)
	351	{
	352	// Analysis from AOD: to be implemented...
	353	// in the near future.
	354	fAOD->Print();
	355
	356	}
	357	//_________________________________________________________________________//
	358	void AliProtonFeedDownAnalysis::Analyze(AliStack *stack)
	359	{
	360	Double_t containerInput[2] ;
	361	Float_t weight;
	362	for (Int_t ipart=0; ipart<stack->GetNtrack(); ipart++)
	363	{
	364	TParticle *particle = stack->Particle(ipart);
	365	Int_t code=particle->GetPdgCode();
c9fcf8c9	366	/*if (code==3122)
198dfda4	367	{
	368	fLambda->Fill(particle->Pt());
	369	fLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));
	370	}
	371	if (code==-3122)
	372	{
	373	fAntiLambda->Fill(particle->Pt());
	374	fAntiLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));
c9fcf8c9	375	}*/
198dfda4	376	if (TMath::Abs(code)==2212)
	377	{
	378	Int_t trackflag=LambdaIsMother(ipart,stack);//1 mother lambda -1 mother anti lambda 0 mother something else
	379	if (trackflag!=0)
	380	weight=GetWeightforProton(ipart,stack);
	381	else
	382	weight=1.0;
	383	if(particle->GetUniqueID() == 13) //recject hadronic inter.
	384	continue;
	385	if((particle->Pt() > fMaxPt)\|\|(particle->Pt() < fMinPt))
	386	continue;
	387	if(fProtonAnalysisBase->GetEtaMode())
	388	{
	389	if((particle->Eta()> fMaxY)\|\|(particle->Eta()< fMinY))
	390	continue;
	391	}
	392	else
	393	{
	394	if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)\|\|(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz())< fMinY))
	395	continue;
	396	}
	397	if(fProtonAnalysisBase->GetEtaMode())
	398	{
	399	containerInput[0] =particle->Eta();
	400	}
	401	else
	402	{
	403	containerInput[0] = fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz());
	404	}
	405	//containerInput[0] = particle->Eta() ;
	406	containerInput[1] = (Float_t)particle->Pt();
	407	if (particle->GetPdgCode()==2212)
	408	{
	409	fProtonContainer->Fill(containerInput,kAll,weight);
	410	if(ipart<stack->GetNprimary())
	411	{
	412	fProtonContainer->Fill(containerInput,kPrimary,weight);
	413	}
	414	else
	415	{
	416	if(trackflag==1)
	417	fProtonContainer->Fill(containerInput,kFromLambda,weight);
	418	}
	419	}
	420	if (particle->GetPdgCode()==-2212)
	421	{
	422	fAntiProtonContainer->Fill(containerInput,kAll,weight);
	423	if(ipart<stack->GetNprimary())
	424	{
	425	fAntiProtonContainer->Fill(containerInput,kPrimary,weight);
	426	}
	427	else
	428	{
	429	if(trackflag==-1)
	430	fAntiProtonContainer->Fill(containerInput,kFromLambda,weight);
	431	}
	432	}
	433	}
	434	}
	435	}
	436	//______________________________________________________________________________________________
	437	Int_t AliProtonFeedDownAnalysis::LambdaIsMother(Int_t number, AliStack *stack)
	438	{
	439	if(number<0)
440	return 0;
441	TParticle *particle = stack->Particle(number);
442	Int_t motherPDG=0;
443	Int_t lmother=-1;
444	lmother=particle->GetFirstMother();
445	if (lmother<0)
446	return 0;
447	TParticle *mparticle=stack->Particle(lmother);
448	motherPDG=mparticle->GetPdgCode();
449	if(motherPDG==3122)
450	return 1;
451	if(motherPDG==-3122)
452	return -1;
453	return 0;
454	}
455	//___________________________________________________________________________________________
456	Float_t AliProtonFeedDownAnalysis::GetWeightforProton(Int_t number,AliStack *stack)
457	{
458	if(number<0)
459	return 1.0;
460	TParticle *particle = stack->Particle(number);
461	Int_t lmother=particle->GetFirstMother();
462	TParticle *mparticle=stack->Particle(lmother);
463	return fweightfunction->Eval(mparticle->Pt());
464	}
465	//______________________________________________________________________________________________
466	Float_t AliProtonFeedDownAnalysis::GetWeightforLambda(Float_t pt)
467	{
468	return fweightfunction->Eval(pt);
469	}