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

//---- AliRoot system ----
#include "AliAnaExample.h"
#include "AliCaloTrackReader.h"
#include "AliAODPWG4Particle.h"
#include "AliESDCaloCells.h"
#include "AliStack.h"
#include "AliCaloPID.h"
#include "AliFiducialCut.h"
#include "AliAODCaloCells.h"
#include "AliAODCaloCluster.h"
#include "AliAODTrack.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
    
  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("PWG4Particle");
  AddToHistogramsName("AnaExample_");

  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("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
  AliAnaPartCorrBaseClass::Print(" ");

  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("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
    if(fDetector == "CTS" && GetAODCTS())printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
    if(fDetector == "PHOS" && GetAODPHOS())printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
  }
  
  //Get List with tracks or clusters  
  TObjArray * partList = new TObjArray();
  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"){
    
    //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 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(IsFiducialCutOn())
	in =  GetFiducialCut()->IsInFiducialCut(mom,fDetector) ;
      
      if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Cluster pt %2.2f, phi %2.2f, pdg %d, in fiducial 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) {
	printf("AliAnaExample::MakeAnalysisFillAOD() - STOP: No CELLS available for analysis");
	abort();
      }
      //Some prints
      if(GetDebug() > 0 && esdCell )
	printf("AliAnaExample::MakeAnalysisFillAOD() - 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("AliAnaExample::MakeAnalysisFillAOD() - 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 =  GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
      if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial 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("AliAnaExample::MakeAnalysisFillAOD() - 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("AliAnaExample::MakeAnalysisFillHistograms() - 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
	  Float_t ptprim  = 0;
	  Float_t phiprim = 0;
	  Float_t etaprim = 0;
	  if(GetReader()->ReadStack()){
		  AliStack * stack =  GetMCStack() ;
      
		  if(ph->GetLabel() < 0 || !stack) {
			  printf("AliAnaExample::MakeAnalysisFillHistograms() *** bad label or no stack ***:  label %d \n", ph->GetLabel());
			  continue;
		  }
      
		  if(ph->GetLabel() >=  stack->GetNtrack()) {
			  printf("AliAnaExample::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
			  continue ;
		  }
      
		  TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
      	  ptprim  = mom->Pt();
		  phiprim = mom->Phi();
		  etaprim = mom->Eta();
	  }
	  else if(GetReader()->ReadAODMCParticles()){
	  	  printf("AliAnaExample::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
	  }
	  fh2Pt->Fill(ph->Pt(), ptprim);
	  fh2Phi->Fill(ph->Phi(), phiprim);
	  fh2Eta->Fill(ph->Eta(), etaprim);
    }//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("AliAnaExample::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
	fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
      }
    }//calo cells container exist
  }
}

//________________________________________________________________________
void AliAnaExample::ReadHistograms(TList* outputList)
{
	// Needed when Terminate is executed in distributed environment
	// Refill analysis histograms of this class with corresponding histograms in output list. 
	
	// Histograms of this analsys are kept in the same list as other analysis, recover the position of
	// the first one and then add the next 
	Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hPt"));
	
	//Read histograms, must be in the same order as in GetCreateOutputObject.
	fhPt   = (TH1F *) outputList->At(index++); 
	fhPhi  = (TH1F *) outputList->At(index++); 
	fhEta  = (TH1F *) outputList->At(index++);
	
	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
	 fhNCells     = (TH1F *) outputList->At(index++); 
	 fhAmplitude  = (TH1F *) outputList->At(index++); 
	}
	
	if(IsDataMC()){
	  fh2Pt  = (TH2F *) outputList->At(index++); 
	  fh2Phi = (TH2F *) outputList->At(index++); 
	  fh2Eta = (TH2F *) outputList->At(index); 
	}
	
}

//__________________________________________________________________
void  AliAnaExample::Terminate(TList* outputList) 
{
  
  //Do some plots to end
  
  //Recover histograms from output histograms list, needed for distributed analysis.	
  ReadHistograms(outputList);
	
  printf(" AliAnaExample::Terminate()  *** %s Report:", GetName()) ; 
  printf(" AliAnaExample::Terminate()        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("AliAnaExample::Terminate() - !! 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"
591cc579	31	#include "TObjArray.h"
1c5acb87	32	#include "TParticle.h"
6639984f	33	#include "TCanvas.h"
	34	#include "TPad.h"
	35	#include "TROOT.h"
477d6cee	36	#include "TH2F.h"
1c5acb87	37
	38	//---- AliRoot system ----
	39	#include "AliAnaExample.h"
	40	#include "AliCaloTrackReader.h"
1c5acb87	41	#include "AliAODPWG4Particle.h"
	42	#include "AliESDCaloCells.h"
	43	#include "AliStack.h"
	44	#include "AliCaloPID.h"
ff45398a	45	#include "AliFiducialCut.h"
477d6cee	46	#include "AliAODCaloCells.h"
	47	#include "AliAODCaloCluster.h"
	48	#include "AliAODTrack.h"
1c5acb87	49
	50	ClassImp(AliAnaExample)
	51
	52	//____________________________________________________________________________
	53	AliAnaExample::AliAnaExample() :
	54	AliAnaPartCorrBaseClass(),fPdg(0), fDetector(""), fhPt(0),fhPhi(0),fhEta(0), fh2Pt(0),fh2Phi(0),fh2Eta(0),
	55	fhNCells(0), fhAmplitude(0)
	56	{
	57	//Default Ctor
	58
	59	//Initialize parameters
	60	InitParameters();
	61	}
	62
	63	//____________________________________________________________________________
	64	AliAnaExample::AliAnaExample(const AliAnaExample & ex) :
	65	AliAnaPartCorrBaseClass(ex), fPdg(ex.fPdg), fDetector(ex.fDetector), fhPt(ex.fhPt), fhPhi(ex.fhPhi),fhEta(ex.fhEta),
	66	fh2Pt(ex.fh2Pt), fh2Phi(ex.fh2Phi),fh2Eta(ex.fh2Eta), fhNCells(ex.fhNCells), fhAmplitude(ex.fhAmplitude)
	67	{
	68	// cpy ctor
	69
	70	}
	71
	72	//_________________________________________________________________________
	73	AliAnaExample & AliAnaExample::operator = (const AliAnaExample & ex)
	74	{
	75	// assignment operator
	76
	77	if(this == &ex)return *this;
	78	((AliAnaPartCorrBaseClass *)this)->operator=(ex);
	79
	80	fPdg = ex.fPdg;
	81	fDetector = ex.fDetector;
	82	fhPt = ex.fhPt;
	83	fhPhi = ex.fhPhi;
	84	fhEta = ex.fhEta;
	85	fh2Pt = ex.fh2Pt;
	86	fh2Phi = ex.fh2Phi;
	87	fh2Eta = ex.fh2Eta;
	88	fhNCells = ex.fhNCells;
	89	fhAmplitude = ex.fhAmplitude;
	90
	91	return *this;
	92
	93	}
	94
	95	// //____________________________________________________________________________
	96	// AliAnaExample::~AliAnaExample()
	97	// {
	98	// // Remove all pointers except analysis output pointers.
	99
	100	// ;
	101	// }
	102
	103
	104	//________________________________________________________________________
	105	TList * AliAnaExample::GetCreateOutputObjects()
	106	{
477d6cee	107	// Create histograms to be saved in output file and
	108	// store them in fOutputContainer
	109
	110	TList * outputContainer = new TList() ;
	111	outputContainer->SetName("ExampleHistos") ;
	112
	113	Int_t nptbins = GetHistoNPtBins();
	114	Int_t nphibins = GetHistoNPhiBins();
	115	Int_t netabins = GetHistoNEtaBins();
	116	Float_t ptmax = GetHistoPtMax();
	117	Float_t phimax = GetHistoPhiMax();
	118	Float_t etamax = GetHistoEtaMax();
	119	Float_t ptmin = GetHistoPtMin();
	120	Float_t phimin = GetHistoPhiMin();
	121	Float_t etamin = GetHistoEtaMin();
	122
	123	fhPt = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
	124	fhPt->SetXTitle("p_{T} (GeV/c)");
	125	outputContainer->Add(fhPt);
	126
	127	fhPhi = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
	128	fhPhi->SetXTitle("#phi (rad)");
	129	outputContainer->Add(fhPhi);
	130
	131	fhEta = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
	132	fhEta->SetXTitle("#eta ");
	133	outputContainer->Add(fhEta);
	134
	135	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
	136
	137	//Calo cells
	138	fhNCells = new TH1F ("hNCells","# cells per event", 100,0,1000);
	139	fhNCells->SetXTitle("n cells");
	140	outputContainer->Add(fhNCells);
	141
	142	fhAmplitude = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
	143	fhAmplitude->SetXTitle("Amplitude ");
	144	outputContainer->Add(fhAmplitude);
	145	}
	146
	147	if(IsDataMC()){
	148	fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
	149	fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
	150	fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
	151	outputContainer->Add(fh2Pt);
	152
	153	fh2Phi = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
	154	fh2Phi->SetXTitle("#phi_{rec} (rad)");
	155	fh2Phi->SetYTitle("#phi_{gen} (rad)");
	156	outputContainer->Add(fh2Phi);
	157
	158	fh2Eta = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
	159	fh2Eta->SetXTitle("#eta_{rec} ");
	160	fh2Eta->SetYTitle("#eta_{gen} ");
	161	outputContainer->Add(fh2Eta);
	162
	163	}
a3aebfff	164
477d6cee	165	return outputContainer;
1c5acb87	166	}
1c5acb87	167
477d6cee	168	//__________________________________________________
1c5acb87	169	void AliAnaExample::InitParameters()
	170	{
	171	//Initialize the parameters of the analysis.
	172	SetOutputAODClassName("AliAODPWG4Particle");
a3aebfff	173	SetOutputAODName("PWG4Particle");
	174	AddToHistogramsName("AnaExample_");
	175
1c5acb87	176	fPdg = 22; //Keep photons
	177	fDetector = "PHOS";
	178
	179	}
	180
	181	//__________________________________________________________________
	182	void AliAnaExample::Print(const Option_t * opt) const
	183	{
	184	//Print some relevant parameters set for the analysis
	185	if(! opt)
	186	return;
	187
a3aebfff	188	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
	189	AliAnaPartCorrBaseClass::Print(" ");
	190
1c5acb87	191	printf("Select clusters with pdg %d \n",fPdg);
	192	printf("Select detector %s \n",fDetector.Data());
	193
	194	}
	195
	196	//__________________________________________________________________
	197	void AliAnaExample::MakeAnalysisFillAOD()
	198	{
477d6cee	199	//Do analysis and fill aods
	200
	201	//Some prints
	202	if(GetDebug() > 0){
a3aebfff	203	if(fDetector == "EMCAL" && GetAODEMCAL())printf("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
	204	if(fDetector == "CTS" && GetAODCTS())printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
	205	if(fDetector == "PHOS" && GetAODPHOS())printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
477d6cee	206	}
	207
	208	//Get List with tracks or clusters
591cc579	209	TObjArray * partList = new TObjArray();
477d6cee	210	if(fDetector == "CTS") partList = GetAODCTS();
	211	else if(fDetector == "EMCAL") partList = GetAODEMCAL();
	212	else if(fDetector == "PHOS") partList = GetAODPHOS();
	213
	214	if(!partList \|\| partList->GetEntriesFast() == 0) return ;
	215
	216	//Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
	217	if(fDetector == "EMCAL" \|\| fDetector == "PHOS"){
	218
	219	//Get vertex for photon momentum calculation
	220	Double_t v[3] ; //vertex ;
	221	GetReader()->GetVertex(v);
	222
	223	TLorentzVector mom ;
	224	for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
	225
	226	AliAODCaloCluster * calo = (AliAODCaloCluster*) (partList->At(i));
	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;
ff45398a	243	if(IsFiducialCutOn())
ff45398a	244	in = GetFiducialCut()->IsInFiducialCut(mom,fDetector) ;
477d6cee	245
ff45398a	246	if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Cluster pt %2.2f, phi %2.2f, pdg %d, in fiducial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
477d6cee	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") {
15800db4	265	ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
477d6cee	266	esdCell = (AliESDCaloCells *) GetPHOSCells();
	267	}
	268	else {
15800db4	269	ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
477d6cee	270	esdCell = (AliESDCaloCells *) GetEMCALCells();
	271	}
	272
	273	if(!esdCell) {
591cc579	274	printf("AliAnaExample::MakeAnalysisFillAOD() - STOP: No CELLS available for analysis");
477d6cee	275	abort();
	276	}
	277	//Some prints
	278	if(GetDebug() > 0 && esdCell )
a3aebfff	279	printf("AliAnaExample::MakeAnalysisFillAOD() - In ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
477d6cee	280
	281	//Fill AODCells in file
	282	Int_t ncells = esdCell->GetNumberOfCells() ;
	283	GetAODCaloCells()->CreateContainer(ncells);
	284
	285	GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
	286
	287	for (Int_t iCell = 0; iCell < ncells; iCell++) {
a3aebfff	288	if(GetDebug() > 2) printf("AliAnaExample::MakeAnalysisFillAOD() - Cell : amp %f, absId %d, time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
1c5acb87	289
477d6cee	290	GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
	291	}
	292	GetAODCaloCells()->Sort();
	293	}
	294	}//cluster-cell analysis
	295	else if(fDetector == "CTS"){ //Track analysis
	296	//Fill AODParticle with CTS aods
	297	TVector3 p3;
	298	for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
	299
	300	AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(i));
	301
	302	//Fill AODParticle after some selection
	303	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
	304	p3.SetXYZ(mom[0],mom[1],mom[2]);
	305
	306	//Acceptance selection
ff45398a	307	Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
ff45398a	308	if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
477d6cee	309	if(p3.Pt() > GetMinPt() && in) {
	310	AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
	311	tr.SetDetector("CTS");
	312	AddAODParticle(tr);
	313	}//selection
	314	}//loop
	315	}//CTS analysis
	316
	317	if(GetDebug() > 0) {
	318	if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
	319	//printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
a3aebfff	320	printf("AliAnaExample::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
477d6cee	321	// if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
	322	//printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
	323	}
1c5acb87	324	}
	325
	326	//__________________________________________________________________
	327	void AliAnaExample::MakeAnalysisFillHistograms()
	328	{
477d6cee	329	//Do analysis and fill histograms
	330
	331	//Loop on stored AODParticles
	332	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff	333	if(GetDebug() > 0) printf("AliAnaExample::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
477d6cee	334	for(Int_t iaod = 0; iaod < naod ; iaod++){
	335	AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
	336
	337	fhPt->Fill(ph->Pt());
	338	fhPhi->Fill(ph->Phi());
	339	fhEta->Fill(ph->Eta());
	340
	341	if(IsDataMC()){
	342	//Play with the MC stack if available
591cc579	343	Float_t ptprim = 0;
	344	Float_t phiprim = 0;
	345	Float_t etaprim = 0;
	346	if(GetReader()->ReadStack()){
	347	AliStack * stack = GetMCStack() ;
477d6cee	348
591cc579	349	if(ph->GetLabel() < 0 \|\| !stack) {
	350	printf("AliAnaExample::MakeAnalysisFillHistograms() * bad label or no stack *: label %d \n", ph->GetLabel());
	351	continue;
	352	}
477d6cee	353
591cc579	354	if(ph->GetLabel() >= stack->GetNtrack()) {
	355	printf("AliAnaExample::MakeAnalysisFillHistograms() * large label *: label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
	356	continue ;
	357	}
477d6cee	358
591cc579	359	TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
	360	ptprim = mom->Pt();
	361	phiprim = mom->Phi();
	362	etaprim = mom->Eta();
	363	}
	364	else if(GetReader()->ReadAODMCParticles()){
	365	printf("AliAnaExample::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
	366	}
	367	fh2Pt->Fill(ph->Pt(), ptprim);
	368	fh2Phi->Fill(ph->Phi(), phiprim);
	369	fh2Eta->Fill(ph->Eta(), etaprim);
477d6cee	370	}//Work with stack also
	371	}// aod branch loop
	372
	373	// CaloCells histograms
	374	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
	375	if(GetAODCaloCells()){
	376
	377	Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
	378	fhNCells->Fill(ncells) ;
	379
	380	for (Int_t iCell = 0; iCell < ncells; iCell++) {
a3aebfff	381	if(GetDebug() > 2) printf("AliAnaExample::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
477d6cee	382	fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
	383	}
	384	}//calo cells container exist
	385	}
1c5acb87	386	}
6639984f	387
a5cc4f03	388	//________________________________________________________________________
	389	void AliAnaExample::ReadHistograms(TList* outputList)
	390	{
	391	// Needed when Terminate is executed in distributed environment
	392	// Refill analysis histograms of this class with corresponding histograms in output list.
	393
	394	// Histograms of this analsys are kept in the same list as other analysis, recover the position of
	395	// the first one and then add the next
	396	Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hPt"));
	397
	398	//Read histograms, must be in the same order as in GetCreateOutputObject.
	399	fhPt = (TH1F *) outputList->At(index++);
	400	fhPhi = (TH1F *) outputList->At(index++);
	401	fhEta = (TH1F *) outputList->At(index++);
	402
	403	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
	404	fhNCells = (TH1F *) outputList->At(index++);
	405	fhAmplitude = (TH1F *) outputList->At(index++);
	406	}
	407
	408	if(IsDataMC()){
	409	fh2Pt = (TH2F *) outputList->At(index++);
	410	fh2Phi = (TH2F *) outputList->At(index++);
	411	fh2Eta = (TH2F *) outputList->At(index);
	412	}
	413
	414	}
	415
6639984f	416	//__________________________________________________________________
a5cc4f03	417	void AliAnaExample::Terminate(TList* outputList)
6639984f	418	{
477d6cee	419
6639984f	420	//Do some plots to end
477d6cee	421
a5cc4f03	422	//Recover histograms from output histograms list, needed for distributed analysis.
	423	ReadHistograms(outputList);
	424
a3aebfff	425	printf(" AliAnaExample::Terminate() *** %s Report:", GetName()) ;
a3aebfff	426	printf(" AliAnaExample::Terminate() pt : %5.3f , RMS : %5.3f \n", fhPt->GetMean(), fhPt->GetRMS() ) ;
6639984f	427
	428	TCanvas * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
	429	c->Divide(1, 3);
	430
	431	c->cd(1) ;
	432	gPad->SetLogy();
	433	fhPt->SetLineColor(2);
	434	fhPt->Draw();
	435
	436	c->cd(2) ;
	437	fhPhi->SetLineColor(2);
	438	fhPhi->Draw();
	439
	440	c->cd(3) ;
	441	fhEta->SetLineColor(2);
	442	fhEta->Draw();
	443
	444	c->Print("Example.eps");
	445
	446	char line[1024] ;
	447	sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
	448	gROOT->ProcessLine(line);
	449	sprintf(line, ".!rm -fR *.eps");
	450	gROOT->ProcessLine(line);
	451
a3aebfff	452	printf("AliAnaExample::Terminate() - !! All the eps files are in %s.tar.gz !!!", GetName());
6639984f	453
6639984f	454	}