[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaExample.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.                  *
 **************************************************************************/
/* $Id: $ */

//_________________________________________________________________________
// Example class on how to read AODCaloClusters, ESDCaloCells and AODTracks and how 
// fill AODs with PWG4PartCorr analysis frame
// Select the type of detector information that you want to analyze, CTS (tracking), PHOS or EMCAL
// Select the PID custer type of the calorimeters
// Set min momentum of the cluster/tracks
// Fill few histograms
//
//-- Author: Gustavo Conesa (INFN-LNF)
//_________________________________________________________________________


// --- ROOT system ---
//#include "Riostream.h"
#include "TClonesArray.h"
#include "TParticle.h"

//---- AliRoot system ----
#include "AliAnaExample.h"
#include "AliCaloTrackReader.h"
#include "AliLog.h"
#include "AliAODPWG4Particle.h"
#include "AliESDCaloCells.h"
#include "AliStack.h"
#include "AliCaloPID.h"

ClassImp(AliAnaExample)
  
//____________________________________________________________________________
  AliAnaExample::AliAnaExample() : 
    AliAnaPartCorrBaseClass(),fPdg(0),  fDetector(""), fhPt(0),fhPhi(0),fhEta(0),  fh2Pt(0),fh2Phi(0),fh2Eta(0),
    fhNCells(0), fhAmplitude(0)
{
  //Default Ctor

  //Initialize parameters
  InitParameters();
}

//____________________________________________________________________________
AliAnaExample::AliAnaExample(const AliAnaExample & ex) :   
  AliAnaPartCorrBaseClass(ex), fPdg(ex.fPdg), fDetector(ex.fDetector), fhPt(ex.fhPt),  fhPhi(ex.fhPhi),fhEta(ex.fhEta), 
  fh2Pt(ex.fh2Pt),  fh2Phi(ex.fh2Phi),fh2Eta(ex.fh2Eta), fhNCells(ex.fhNCells), fhAmplitude(ex.fhAmplitude)
{
  // cpy ctor
  
}

//_________________________________________________________________________
AliAnaExample & AliAnaExample::operator = (const AliAnaExample & ex)
{
  // assignment operator

  if(this == &ex)return *this;
  ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
 
  fPdg = ex.fPdg;
  fDetector = ex.fDetector;
  fhPt = ex.fhPt;
  fhPhi = ex.fhPhi;
  fhEta = ex.fhEta;
  fh2Pt = ex.fh2Pt;
  fh2Phi = ex.fh2Phi;
  fh2Eta = ex.fh2Eta;
  fhNCells = ex.fhNCells;
  fhAmplitude = ex.fhAmplitude;

  return *this;

}

// //____________________________________________________________________________
// AliAnaExample::~AliAnaExample() 
// {
//   // Remove all pointers except analysis output pointers.

//   ;
// }


//________________________________________________________________________
TList *  AliAnaExample::GetCreateOutputObjects()
{  
	// Create histograms to be saved in output file and 
	// store them in fOutputContainer
	
	AliDebug(1,"Init parton histograms");
	
	TList * outputContainer = new TList() ; 
	outputContainer->SetName("ExampleHistos") ; 
	
	Int_t nptbins  = GetHistoNPtBins();
	Int_t nphibins = GetHistoNPhiBins();
	Int_t netabins = GetHistoNEtaBins();
	Float_t ptmax  = GetHistoPtMax();
	Float_t phimax = GetHistoPhiMax();
	Float_t etamax = GetHistoEtaMax();
	Float_t ptmin  = GetHistoPtMin();
	Float_t phimin = GetHistoPhiMin();
	Float_t etamin = GetHistoEtaMin();	
	
	fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax); 
	fhPt->SetXTitle("p_{T} (GeV/c)");
	outputContainer->Add(fhPt);
	
	fhPhi  = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax); 
	fhPhi->SetXTitle("#phi (rad)");
	outputContainer->Add(fhPhi);
	
	fhEta  = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax); 
	fhEta->SetXTitle("#eta ");
	outputContainer->Add(fhEta);
	
	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
		
		//Calo cells
		fhNCells  = new TH1F ("hNCells","# cells per event", 100,0,1000); 
		fhNCells->SetXTitle("n cells");
		outputContainer->Add(fhNCells);
		
		fhAmplitude  = new TH1F ("hAmplitude","#eta distribution", 100,0,1000); 
		fhAmplitude->SetXTitle("Amplitude ");
		outputContainer->Add(fhAmplitude);
	} 
	
	if(IsDataMC()){
		fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
		fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
		fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
		outputContainer->Add(fh2Pt);
		
		fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax); 
		fh2Phi->SetXTitle("#phi_{rec} (rad)");
		fh2Phi->SetYTitle("#phi_{gen} (rad)");
		outputContainer->Add(fh2Phi);
		
		fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax); 
		fh2Eta->SetXTitle("#eta_{rec} ");
		fh2Eta->SetYTitle("#eta_{gen} ");
		outputContainer->Add(fh2Eta);
		
	}
	return outputContainer;
}

 //__________________________________________________
void AliAnaExample::InitParameters()
{ 
  //Initialize the parameters of the analysis.
  SetOutputAODClassName("AliAODPWG4Particle");
  SetOutputAODName("example");
  fPdg = 22; //Keep photons
  fDetector = "PHOS";
  
}

//__________________________________________________________________
void AliAnaExample::Print(const Option_t * opt) const
{
  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  printf("Min Pt = %3.2f\n", GetMinPt());
  printf("Max Pt = %3.2f\n", GetMaxPt());
  printf("Select clusters with pdg %d \n",fPdg);
  printf("Select detector %s \n",fDetector.Data());
  
} 

//__________________________________________________________________
void  AliAnaExample::MakeAnalysisFillAOD() 
{
	//Do analysis and fill aods
	
	//Some prints
	if(GetDebug() > 0){
		if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
		if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
		if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
	}
	
	//Get List with tracks or clusters  
	TClonesArray * partList = new TClonesArray;
	if(fDetector == "CTS") partList = GetAODCTS();
	else if(fDetector == "EMCAL") partList = GetAODEMCAL();
	else if(fDetector == "PHOS") partList = GetAODPHOS();
	
	if(!partList || partList->GetEntriesFast() == 0) return ;
	
	//Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
	if(fDetector == "EMCAL" || fDetector == "PHOS"){
		
		//WORK WITH CALOCLUSTERS 
		if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC){ 
			ConnectAODCaloClusters(); //Do Only when filling AODCaloClusters 
			if(GetDebug() > 0) printf("Example: in calo clusters aod entries %d\n", GetAODCaloClusters()->GetEntriesFast());
		}
		//Get vertex for photon momentum calculation
		Double_t v[3] ; //vertex ;
		GetReader()->GetVertex(v);
		
		TLorentzVector mom ;
		for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
			
			AliAODCaloCluster * calo =  (AliAODCaloCluster*) (partList->At(i));
			
			//Fill AODCaloClusters  
			if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
				AddAODCaloCluster(AliAODCaloCluster(*(calo)));
			
			//Fill AODParticle after some selection
			calo->GetMomentum(mom,v);
			Int_t pdg = fPdg;
			
			if(IsCaloPIDOn()){
				Double_t pid[13];
				calo->GetPID(pid);
				pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
				//pdg = GetCaloPID()->GetPdg(fDetector,mom,
				//		  calo->GetM02(), calo->GetM02(),
				//		  calo->GetDispersion(), 0, 0); 
			}
			
			//Acceptance selection   
			Bool_t in = kTRUE;
			if(IsFidutialCutOn())
				in =  GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
			
			if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
			
			//Select cluster if momentum, pdg and acceptance are good
			if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
				AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
				//AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
				ph.SetLabel(calo->GetLabel(0));
				ph.SetPdg(pdg);
				ph.SetDetector(fDetector);
				AddAODParticle(ph);
			}//selection
		}//loop
		
		if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
			//WORK WITH ESDCALOCELLS
			//Don't connect in the same analysis PHOS and EMCAL cells.
			
			AliESDCaloCells * esdCell = new AliESDCaloCells ;
			if(fDetector == "PHOS") {
				ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
				esdCell = (AliESDCaloCells *) GetPHOSCells();
			}
			else  {
				ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
				esdCell = (AliESDCaloCells *) GetEMCALCells();
			}
			
			if(!esdCell) AliFatal("No CELLS available for analysis");
			
			//Some prints
			if(GetDebug() > 0 && esdCell )
				printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());    
			
			//Fill AODCells in file
			Int_t ncells = esdCell->GetNumberOfCells() ;
			GetAODCaloCells()->CreateContainer(ncells);
			
			GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
			
			for (Int_t iCell = 0; iCell < ncells; iCell++) {      
				if(GetDebug() > 2)  printf("cell : amp %f, absId %d,  time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
				
				GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
			}
			GetAODCaloCells()->Sort();
		} 
	}//cluster-cell analysis
	else if(fDetector == "CTS"){ //Track analysis
		//Fill AODParticle with CTS aods
		TVector3 p3;
		for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
			
			AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(i));
			
			//Fill AODParticle after some selection       
			Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
			p3.SetXYZ(mom[0],mom[1],mom[2]);
			
			//Acceptance selection
			Bool_t in =  GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
			if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
			if(p3.Pt() > GetMinPt() && in) {
				AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
				tr.SetDetector("CTS");
				AddAODParticle(tr);
			}//selection
		}//loop
	}//CTS analysis
	
	if(GetDebug() > 0) {
		if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
			printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
		printf("Example: final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());  
		if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
			printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
	}
} 

//__________________________________________________________________
void  AliAnaExample::MakeAnalysisFillHistograms() 
{
	//Do analysis and fill histograms
	
	//Loop on stored AODParticles
	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
	if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
	for(Int_t iaod = 0; iaod < naod ; iaod++){
		AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
		
		fhPt->Fill(ph->Pt());
		fhPhi->Fill(ph->Phi());
		fhEta->Fill(ph->Eta());
		
		if(IsDataMC()){
			//Play with the MC stack if available
			AliStack * stack =  GetMCStack() ;
			
			if(ph->GetLabel() < 0 || !stack) {
				printf("*** bad label or no stack ***:  label %d \n", ph->GetLabel());
				continue;
			}
			
			if(ph->GetLabel() >=  stack->GetNtrack()) {
				printf("*** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
				continue ;
			}
			
			TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
			
			fh2Pt->Fill(ph->Pt(), mom->Pt());
			fh2Phi->Fill(ph->Phi(), mom->Phi());
			fh2Eta->Fill(ph->Eta(), mom->Eta());
		}//Work with stack also
	}// aod branch loop
	
	// CaloCells histograms
	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
		if(GetAODCaloCells()){
			
			Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
			fhNCells->Fill(ncells) ;
			
			for (Int_t iCell = 0; iCell < ncells; iCell++) {      
				if(GetDebug() > 2)  printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
				fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
			}
		}//calo cells container exist
	}
}
Commit	Line	Data
1c5acb87	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	/* $Id: $ */
	16
	17	//_________________________________________________________________________
	18	// Example class on how to read AODCaloClusters, ESDCaloCells and AODTracks and how
	19	// fill AODs with PWG4PartCorr analysis frame
	20	// Select the type of detector information that you want to analyze, CTS (tracking), PHOS or EMCAL
	21	// Select the PID custer type of the calorimeters
	22	// Set min momentum of the cluster/tracks
	23	// Fill few histograms
	24	//
	25	//-- Author: Gustavo Conesa (INFN-LNF)
	26	//_________________________________________________________________________
	27
	28
	29	// --- ROOT system ---
	30	//#include "Riostream.h"
	31	#include "TClonesArray.h"
	32	#include "TParticle.h"
	33
	34	//---- AliRoot system ----
	35	#include "AliAnaExample.h"
	36	#include "AliCaloTrackReader.h"
	37	#include "AliLog.h"
	38	#include "AliAODPWG4Particle.h"
	39	#include "AliESDCaloCells.h"
	40	#include "AliStack.h"
	41	#include "AliCaloPID.h"
	42
	43	ClassImp(AliAnaExample)
	44
	45	//____________________________________________________________________________
	46	AliAnaExample::AliAnaExample() :
	47	AliAnaPartCorrBaseClass(),fPdg(0), fDetector(""), fhPt(0),fhPhi(0),fhEta(0), fh2Pt(0),fh2Phi(0),fh2Eta(0),
	48	fhNCells(0), fhAmplitude(0)
	49	{
	50	//Default Ctor
	51
	52	//Initialize parameters
	53	InitParameters();
	54	}
	55
	56	//____________________________________________________________________________
	57	AliAnaExample::AliAnaExample(const AliAnaExample & ex) :
	58	AliAnaPartCorrBaseClass(ex), fPdg(ex.fPdg), fDetector(ex.fDetector), fhPt(ex.fhPt), fhPhi(ex.fhPhi),fhEta(ex.fhEta),
	59	fh2Pt(ex.fh2Pt), fh2Phi(ex.fh2Phi),fh2Eta(ex.fh2Eta), fhNCells(ex.fhNCells), fhAmplitude(ex.fhAmplitude)
	60	{
	61	// cpy ctor
	62
	63	}
	64
65	//_________________________________________________________________________
66	AliAnaExample & AliAnaExample::operator = (const AliAnaExample & ex)
67	{
68	// assignment operator
69
70	if(this == &ex)return *this;
71	((AliAnaPartCorrBaseClass *)this)->operator=(ex);
72
73	fPdg = ex.fPdg;
74	fDetector = ex.fDetector;
75	fhPt = ex.fhPt;
76	fhPhi = ex.fhPhi;
77	fhEta = ex.fhEta;
78	fh2Pt = ex.fh2Pt;
79	fh2Phi = ex.fh2Phi;
80	fh2Eta = ex.fh2Eta;
81	fhNCells = ex.fhNCells;
82	fhAmplitude = ex.fhAmplitude;
83
84	return *this;
85
86	}
87
88	// //____________________________________________________________________________
89	// AliAnaExample::~AliAnaExample()
90	// {
91	// // Remove all pointers except analysis output pointers.
92
93	// ;
94	// }
95
96
97	//________________________________________________________________________
98	TList * AliAnaExample::GetCreateOutputObjects()
99	{
100	// Create histograms to be saved in output file and
101	// store them in fOutputContainer
102
103	AliDebug(1,"Init parton histograms");
104
105	TList * outputContainer = new TList() ;
106	outputContainer->SetName("ExampleHistos") ;
107
108	Int_t nptbins = GetHistoNPtBins();
109	Int_t nphibins = GetHistoNPhiBins();
110	Int_t netabins = GetHistoNEtaBins();
111	Float_t ptmax = GetHistoPtMax();
112	Float_t phimax = GetHistoPhiMax();
113	Float_t etamax = GetHistoEtaMax();
114	Float_t ptmin = GetHistoPtMin();
115	Float_t phimin = GetHistoPhiMin();
116	Float_t etamin = GetHistoEtaMin();
117
118	fhPt = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
119	fhPt->SetXTitle("p_{T} (GeV/c)");
120	outputContainer->Add(fhPt);
121
122	fhPhi = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
123	fhPhi->SetXTitle("#phi (rad)");
124	outputContainer->Add(fhPhi);
125
126	fhEta = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
127	fhEta->SetXTitle("#eta ");
128	outputContainer->Add(fhEta);
129
130	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
131
132	//Calo cells
133	fhNCells = new TH1F ("hNCells","# cells per event", 100,0,1000);
134	fhNCells->SetXTitle("n cells");
135	outputContainer->Add(fhNCells);
136
137	fhAmplitude = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
138	fhAmplitude->SetXTitle("Amplitude ");
139	outputContainer->Add(fhAmplitude);
140	}
141
142	if(IsDataMC()){
143	fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
144	fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
145	fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
146	outputContainer->Add(fh2Pt);
147
148	fh2Phi = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
149	fh2Phi->SetXTitle("#phi_{rec} (rad)");
150	fh2Phi->SetYTitle("#phi_{gen} (rad)");
151	outputContainer->Add(fh2Phi);
152
153	fh2Eta = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
154	fh2Eta->SetXTitle("#eta_{rec} ");
155	fh2Eta->SetYTitle("#eta_{gen} ");
156	outputContainer->Add(fh2Eta);
157
158	}
159	return outputContainer;
160	}
161
162	//__________________________________________________
163	void AliAnaExample::InitParameters()
164	{
165	//Initialize the parameters of the analysis.
166	SetOutputAODClassName("AliAODPWG4Particle");
167	SetOutputAODName("example");
168	fPdg = 22; //Keep photons
169	fDetector = "PHOS";
170
171	}
172
173	//__________________________________________________________________
174	void AliAnaExample::Print(const Option_t * opt) const
175	{
176	//Print some relevant parameters set for the analysis
177	if(! opt)
178	return;
179
180	printf("Min Pt = %3.2f\n", GetMinPt());
181	printf("Max Pt = %3.2f\n", GetMaxPt());
182	printf("Select clusters with pdg %d \n",fPdg);
183	printf("Select detector %s \n",fDetector.Data());
184
185	}
186
187	//__________________________________________________________________
188	void AliAnaExample::MakeAnalysisFillAOD()
189	{
190	//Do analysis and fill aods
191
192	//Some prints
193	if(GetDebug() > 0){
194	if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
195	if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
196	if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
197	}
198
199	//Get List with tracks or clusters
200	TClonesArray * partList = new TClonesArray;
201	if(fDetector == "CTS") partList = GetAODCTS();
202	else if(fDetector == "EMCAL") partList = GetAODEMCAL();
203	else if(fDetector == "PHOS") partList = GetAODPHOS();
204
205	if(!partList \|\| partList->GetEntriesFast() == 0) return ;
206
207	//Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
208	if(fDetector == "EMCAL" \|\| fDetector == "PHOS"){
209
210	//WORK WITH CALOCLUSTERS
211	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC){
212	ConnectAODCaloClusters(); //Do Only when filling AODCaloClusters
213	if(GetDebug() > 0) printf("Example: in calo clusters aod entries %d\n", GetAODCaloClusters()->GetEntriesFast());
214	}
215	//Get vertex for photon momentum calculation
216	Double_t v[3] ; //vertex ;
217	GetReader()->GetVertex(v);
218
219	TLorentzVector mom ;
220	for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
221
222	AliAODCaloCluster * calo = (AliAODCaloCluster*) (partList->At(i));
223
224	//Fill AODCaloClusters
225	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
226	AddAODCaloCluster(AliAODCaloCluster(*(calo)));
227
228	//Fill AODParticle after some selection
229	calo->GetMomentum(mom,v);
230	Int_t pdg = fPdg;
231
232	if(IsCaloPIDOn()){
233	Double_t pid[13];
234	calo->GetPID(pid);
235	pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
236	//pdg = GetCaloPID()->GetPdg(fDetector,mom,
237	// calo->GetM02(), calo->GetM02(),
238	// calo->GetDispersion(), 0, 0);
239	}
240
241	//Acceptance selection
242	Bool_t in = kTRUE;
243	if(IsFidutialCutOn())
244	in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
245
246	if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
247
248	//Select cluster if momentum, pdg and acceptance are good
249	if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
250	AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
251	//AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
252	ph.SetLabel(calo->GetLabel(0));
253	ph.SetPdg(pdg);
254	ph.SetDetector(fDetector);
255	AddAODParticle(ph);
256	}//selection
257	}//loop
258
259	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
260	//WORK WITH ESDCALOCELLS
261	//Don't connect in the same analysis PHOS and EMCAL cells.
262
263	AliESDCaloCells * esdCell = new AliESDCaloCells ;
264	if(fDetector == "PHOS") {
265	ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
266	esdCell = (AliESDCaloCells *) GetPHOSCells();
267	}
268	else {
269	ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
270	esdCell = (AliESDCaloCells *) GetEMCALCells();
271	}
272
273	if(!esdCell) AliFatal("No CELLS available for analysis");
274
275	//Some prints
276	if(GetDebug() > 0 && esdCell )
277	printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
278
279	//Fill AODCells in file
280	Int_t ncells = esdCell->GetNumberOfCells() ;
281	GetAODCaloCells()->CreateContainer(ncells);
282
283	GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
284
285	for (Int_t iCell = 0; iCell < ncells; iCell++) {
286	if(GetDebug() > 2) printf("cell : amp %f, absId %d, time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
287
288	GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
289	}
290	GetAODCaloCells()->Sort();
291	}
292	}//cluster-cell analysis
293	else if(fDetector == "CTS"){ //Track analysis
294	//Fill AODParticle with CTS aods
295	TVector3 p3;
296	for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
297
298	AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(i));
299
300	//Fill AODParticle after some selection
301	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
302	p3.SetXYZ(mom[0],mom[1],mom[2]);
303
304	//Acceptance selection
305	Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
306	if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
307	if(p3.Pt() > GetMinPt() && in) {
308	AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
309	tr.SetDetector("CTS");
310	AddAODParticle(tr);
311	}//selection
312	}//loop
313	}//CTS analysis
314
315	if(GetDebug() > 0) {
316	if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
317	printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
318	printf("Example: final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
319	if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
320	printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
321	}
322	}
323
324	//__________________________________________________________________
325	void AliAnaExample::MakeAnalysisFillHistograms()
326	{
327	//Do analysis and fill histograms
328
329	//Loop on stored AODParticles
330	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
331	if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
332	for(Int_t iaod = 0; iaod < naod ; iaod++){
333	AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
334
335	fhPt->Fill(ph->Pt());
336	fhPhi->Fill(ph->Phi());
337	fhEta->Fill(ph->Eta());
338
339	if(IsDataMC()){
340	//Play with the MC stack if available
341	AliStack * stack = GetMCStack() ;
342
343	if(ph->GetLabel() < 0 \|\| !stack) {
344	printf("* bad label or no stack *: label %d \n", ph->GetLabel());
345	continue;
346	}
347
348	if(ph->GetLabel() >= stack->GetNtrack()) {
349	printf("* large label *: label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
350	continue ;
351	}
352
353	TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
354
355	fh2Pt->Fill(ph->Pt(), mom->Pt());
356	fh2Phi->Fill(ph->Phi(), mom->Phi());
357	fh2Eta->Fill(ph->Eta(), mom->Eta());
358	}//Work with stack also
359	}// aod branch loop
360
361	// CaloCells histograms
362	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
363	if(GetAODCaloCells()){
364
365	Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
366	fhNCells->Fill(ncells) ;
367
368	for (Int_t iCell = 0; iCell < ncells; iCell++) {
369	if(GetDebug() > 2) printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
370	fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
371	}
372	}//calo cells container exist
373	}
374	}