[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"
#include "TCanvas.h"
#include "TPad.h"
#include "TROOT.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
	}
}

//__________________________________________________________________
void  AliAnaExample::Terminate() 
{

  //Do some plots to end

  printf(" *** %s Report:", GetName()) ; 
  printf("        pt         : %5.3f , RMS : %5.3f \n", fhPt->GetMean(),   fhPt->GetRMS() ) ;
 
  TCanvas  * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
  c->Divide(1, 3);

  c->cd(1) ; 
  gPad->SetLogy();
  fhPt->SetLineColor(2);
  fhPt->Draw();

  c->cd(2) ; 
  fhPhi->SetLineColor(2);
  fhPhi->Draw();

  c->cd(3) ; 
  fhEta->SetLineColor(2);
  fhEta->Draw();
 
  c->Print("Example.eps");
 
  char line[1024] ; 
  sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ; 
  gROOT->ProcessLine(line);
  sprintf(line, ".!rm -fR *.eps"); 
  gROOT->ProcessLine(line);
 
  printf("!! All the eps files are in %s.tar.gz !!!", GetName());

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