[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 CaloClusters, 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 "AliStack.h"
#include "AliCaloPID.h"
#include "AliFiducialCut.h"
#include "AliVCluster.h"
#include "AliVTrack.h"

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

  //Initialize parameters
  InitParameters();
}

//________________________________________________________________________
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  = GetHistoPtBins();
  Int_t nphibins = GetHistoPhiBins();
  Int_t netabins = GetHistoEtaBins();
  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(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" && GetEMCALClusters())printf("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
    if(fDetector == "CTS"   && GetCTSTracks())  printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n",   GetCTSTracks()  ->GetEntriesFast());
    if(fDetector == "PHOS"  && GetPHOSClusters()) printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n",  GetPHOSClusters() ->GetEntriesFast());
  }
  
  //Get List with tracks or clusters  
  TObjArray * partList = 0x0;
  if(fDetector == "CTS")        partList = GetCTSTracks();
  else if(fDetector == "EMCAL") partList = GetEMCALClusters();
  else if(fDetector == "PHOS")  partList = GetPHOSClusters();
  
  if(!partList || partList->GetEntriesFast() == 0) return ;
  
  //Fill AODParticle with PHOS/EMCAL aods
  if(fDetector == "EMCAL" || fDetector == "PHOS"){
    
    //Get vertex for photon momentum calculation

    
    TLorentzVector mom ;
    for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
      
      AliVCluster * calo =  (AliVCluster*) (partList->At(i));
      
      //Get the index where the cluster comes, to retrieve the corresponding vertex
      Int_t evtIndex = 0 ; 
      if (GetMixedEvent()) {
        evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
      }

      //Fill AODParticle after some selection
      if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
        calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
      else{
        Double_t vertex[]={0,0,0};
        calo->GetMomentum(mom,vertex) ;
      }
      
      Int_t pdg = fPdg;
      
      if(IsCaloPIDOn()){
        const Double_t *pid = calo->GetPID();
        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());
        ph.SetPdg(pdg);
        ph.SetDetector(fDetector);
        AddAODParticle(ph);
      }//selection
    }//loop
  }//Calorimeters
  else if(fDetector == "CTS"){ //Track analysis
    //Fill AODParticle with CTS aods
    TVector3 p3;
    for(Int_t i = 0; i < GetCTSTracks()->GetEntriesFast(); i++){
      
      AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->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("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
  
 }

//________________________________________________________________________
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(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");
  
  const Int_t buffersize = 1024;
  char line[buffersize] ; 
  snprintf(line,buffersize, ".!tar -zcf %s.tar.gz *.eps", GetName()) ; 
  gROOT->ProcessLine(line);
  snprintf(line,buffersize, ".!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	//_________________________________________________________________________
f37fa8d2	18	// Example class on how to read CaloClusters, and AODTracks and how
1c5acb87	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"
1c5acb87	42	#include "AliStack.h"
1c5acb87	43	#include "AliCaloPID.h"
ff45398a	44	#include "AliFiducialCut.h"
f37fa8d2	45	#include "AliVCluster.h"
88f9563f	46	#include "AliVTrack.h"
1c5acb87	47
	48	ClassImp(AliAnaExample)
	49
	50	//____________________________________________________________________________
	51	AliAnaExample::AliAnaExample() :
f37fa8d2	52	AliAnaPartCorrBaseClass(),fPdg(0), fDetector(""), fhPt(0),fhPhi(0),fhEta(0), fh2Pt(0),fh2Phi(0),fh2Eta(0)
1c5acb87	53	{
	54	//Default Ctor
	55
	56	//Initialize parameters
	57	InitParameters();
	58	}
1c5acb87	59
	60	//________________________________________________________________________
	61	TList * AliAnaExample::GetCreateOutputObjects()
	62	{
477d6cee	63	// Create histograms to be saved in output file and
	64	// store them in fOutputContainer
	65
	66	TList * outputContainer = new TList() ;
	67	outputContainer->SetName("ExampleHistos") ;
	68
5a2dbc3c	69	Int_t nptbins = GetHistoPtBins();
	70	Int_t nphibins = GetHistoPhiBins();
	71	Int_t netabins = GetHistoEtaBins();
477d6cee	72	Float_t ptmax = GetHistoPtMax();
	73	Float_t phimax = GetHistoPhiMax();
	74	Float_t etamax = GetHistoEtaMax();
	75	Float_t ptmin = GetHistoPtMin();
	76	Float_t phimin = GetHistoPhiMin();
	77	Float_t etamin = GetHistoEtaMin();
	78
	79	fhPt = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
	80	fhPt->SetXTitle("p_{T} (GeV/c)");
	81	outputContainer->Add(fhPt);
	82
	83	fhPhi = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
	84	fhPhi->SetXTitle("#phi (rad)");
	85	outputContainer->Add(fhPhi);
	86
	87	fhEta = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
	88	fhEta->SetXTitle("#eta ");
	89	outputContainer->Add(fhEta);
	90
477d6cee	91	if(IsDataMC()){
	92	fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
	93	fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
	94	fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
	95	outputContainer->Add(fh2Pt);
	96
	97	fh2Phi = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
	98	fh2Phi->SetXTitle("#phi_{rec} (rad)");
	99	fh2Phi->SetYTitle("#phi_{gen} (rad)");
	100	outputContainer->Add(fh2Phi);
	101
	102	fh2Eta = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
	103	fh2Eta->SetXTitle("#eta_{rec} ");
	104	fh2Eta->SetYTitle("#eta_{gen} ");
	105	outputContainer->Add(fh2Eta);
	106
	107	}
a3aebfff	108
477d6cee	109	return outputContainer;
1c5acb87	110	}
1c5acb87	111
477d6cee	112	//__________________________________________________
1c5acb87	113	void AliAnaExample::InitParameters()
	114	{
	115	//Initialize the parameters of the analysis.
	116	SetOutputAODClassName("AliAODPWG4Particle");
a3aebfff	117	SetOutputAODName("PWG4Particle");
	118	AddToHistogramsName("AnaExample_");
	119
1c5acb87	120	fPdg = 22; //Keep photons
	121	fDetector = "PHOS";
	122
	123	}
	124
	125	//__________________________________________________________________
	126	void AliAnaExample::Print(const Option_t * opt) const
	127	{
	128	//Print some relevant parameters set for the analysis
	129	if(! opt)
	130	return;
	131
a3aebfff	132	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
	133	AliAnaPartCorrBaseClass::Print(" ");
	134
1c5acb87	135	printf("Select clusters with pdg %d \n",fPdg);
	136	printf("Select detector %s \n",fDetector.Data());
	137
	138	}
	139
	140	//__________________________________________________________________
	141	void AliAnaExample::MakeAnalysisFillAOD()
	142	{
477d6cee	143	//Do analysis and fill aods
	144
	145	//Some prints
	146	if(GetDebug() > 0){
be518ab0	147	if(fDetector == "EMCAL" && GetEMCALClusters())printf("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
	148	if(fDetector == "CTS" && GetCTSTracks()) printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks() ->GetEntriesFast());
	149	if(fDetector == "PHOS" && GetPHOSClusters()) printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetPHOSClusters() ->GetEntriesFast());
477d6cee	150	}
	151
	152	//Get List with tracks or clusters
4a745797	153	TObjArray * partList = 0x0;
be518ab0	154	if(fDetector == "CTS") partList = GetCTSTracks();
	155	else if(fDetector == "EMCAL") partList = GetEMCALClusters();
	156	else if(fDetector == "PHOS") partList = GetPHOSClusters();
477d6cee	157
	158	if(!partList \|\| partList->GetEntriesFast() == 0) return ;
	159
0ae57829	160	//Fill AODParticle with PHOS/EMCAL aods
477d6cee	161	if(fDetector == "EMCAL" \|\| fDetector == "PHOS"){
	162
	163	//Get vertex for photon momentum calculation
f8006433	164
477d6cee	165
	166	TLorentzVector mom ;
	167	for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
	168
f37fa8d2	169	AliVCluster * calo = (AliVCluster*) (partList->At(i));
477d6cee	170
f8006433	171	//Get the index where the cluster comes, to retrieve the corresponding vertex
	172	Int_t evtIndex = 0 ;
	173	if (GetMixedEvent()) {
	174	evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
	175	}
	176
477d6cee	177	//Fill AODParticle after some selection
f8006433	178	if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
	179	calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
	180	else{
	181	Double_t vertex[]={0,0,0};
	182	calo->GetMomentum(mom,vertex) ;
	183	}
	184
477d6cee	185	Int_t pdg = fPdg;
	186
	187	if(IsCaloPIDOn()){
f37fa8d2	188	const Double_t *pid = calo->GetPID();
	189	pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
	190	//pdg = GetCaloPID()->GetPdg(fDetector,mom,
	191	// calo->GetM02(), calo->GetM02(),
	192	// calo->GetDispersion(), 0, 0);
477d6cee	193	}
	194
	195	//Acceptance selection
	196	Bool_t in = kTRUE;
ff45398a	197	if(IsFiducialCutOn())
f37fa8d2	198	in = GetFiducialCut()->IsInFiducialCut(mom,fDetector) ;
477d6cee	199
ff45398a	200	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	201
	202	//Select cluster if momentum, pdg and acceptance are good
	203	if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
f37fa8d2	204	AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
	205	//AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
	206	ph.SetLabel(calo->GetLabel());
	207	ph.SetPdg(pdg);
	208	ph.SetDetector(fDetector);
	209	AddAODParticle(ph);
477d6cee	210	}//selection
477d6cee	211	}//loop
f37fa8d2	212	}//Calorimeters
477d6cee	213	else if(fDetector == "CTS"){ //Track analysis
	214	//Fill AODParticle with CTS aods
	215	TVector3 p3;
be518ab0	216	for(Int_t i = 0; i < GetCTSTracks()->GetEntriesFast(); i++){
477d6cee	217
88f9563f	218	AliVTrack * track = (AliVTrack*) (GetCTSTracks()->At(i));
477d6cee	219
	220	//Fill AODParticle after some selection
	221	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
	222	p3.SetXYZ(mom[0],mom[1],mom[2]);
	223
	224	//Acceptance selection
ff45398a	225	Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
ff45398a	226	if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
477d6cee	227	if(p3.Pt() > GetMinPt() && in) {
f37fa8d2	228	AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
	229	tr.SetDetector("CTS");
	230	AddAODParticle(tr);
477d6cee	231	}//selection
	232	}//loop
	233	}//CTS analysis
	234
	235	if(GetDebug() > 0) {
	236	if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
a3aebfff	237	printf("AliAnaExample::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
477d6cee	238	// if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
	239	//printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
	240	}
1c5acb87	241	}
	242
	243	//__________________________________________________________________
	244	void AliAnaExample::MakeAnalysisFillHistograms()
	245	{
477d6cee	246	//Do analysis and fill histograms
	247
	248	//Loop on stored AODParticles
	249	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff	250	if(GetDebug() > 0) printf("AliAnaExample::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
477d6cee	251	for(Int_t iaod = 0; iaod < naod ; iaod++){
	252	AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
	253
	254	fhPt->Fill(ph->Pt());
	255	fhPhi->Fill(ph->Phi());
	256	fhEta->Fill(ph->Eta());
	257
	258	if(IsDataMC()){
f37fa8d2	259	//Play with the MC stack if available
591cc579	260	Float_t ptprim = 0;
	261	Float_t phiprim = 0;
	262	Float_t etaprim = 0;
	263	if(GetReader()->ReadStack()){
	264	AliStack * stack = GetMCStack() ;
477d6cee	265
591cc579	266	if(ph->GetLabel() < 0 \|\| !stack) {
	267	printf("AliAnaExample::MakeAnalysisFillHistograms() * bad label or no stack *: label %d \n", ph->GetLabel());
	268	continue;
	269	}
477d6cee	270
591cc579	271	if(ph->GetLabel() >= stack->GetNtrack()) {
	272	printf("AliAnaExample::MakeAnalysisFillHistograms() * large label *: label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
	273	continue ;
	274	}
477d6cee	275
591cc579	276	TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
f37fa8d2	277	ptprim = mom->Pt();
591cc579	278	phiprim = mom->Phi();
	279	etaprim = mom->Eta();
	280	}
	281	else if(GetReader()->ReadAODMCParticles()){
	282	printf("AliAnaExample::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
	283	}
	284	fh2Pt->Fill(ph->Pt(), ptprim);
	285	fh2Phi->Fill(ph->Phi(), phiprim);
	286	fh2Eta->Fill(ph->Eta(), etaprim);
477d6cee	287	}//Work with stack also
	288	}// aod branch loop
	289
f37fa8d2	290	}
6639984f	291
a5cc4f03	292	//________________________________________________________________________
	293	void AliAnaExample::ReadHistograms(TList* outputList)
	294	{
	295	// Needed when Terminate is executed in distributed environment
	296	// Refill analysis histograms of this class with corresponding histograms in output list.
	297
	298	// Histograms of this analsys are kept in the same list as other analysis, recover the position of
	299	// the first one and then add the next
	300	Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hPt"));
	301
	302	//Read histograms, must be in the same order as in GetCreateOutputObject.
	303	fhPt = (TH1F *) outputList->At(index++);
	304	fhPhi = (TH1F *) outputList->At(index++);
	305	fhEta = (TH1F *) outputList->At(index++);
	306
a5cc4f03	307	if(IsDataMC()){
	308	fh2Pt = (TH2F *) outputList->At(index++);
	309	fh2Phi = (TH2F *) outputList->At(index++);
	310	fh2Eta = (TH2F *) outputList->At(index);
	311	}
	312
	313	}
	314
6639984f	315	//__________________________________________________________________
a5cc4f03	316	void AliAnaExample::Terminate(TList* outputList)
6639984f	317	{
477d6cee	318
6639984f	319	//Do some plots to end
477d6cee	320
a5cc4f03	321	//Recover histograms from output histograms list, needed for distributed analysis.
	322	ReadHistograms(outputList);
	323
a3aebfff	324	printf(" AliAnaExample::Terminate() *** %s Report:", GetName()) ;
a3aebfff	325	printf(" AliAnaExample::Terminate() pt : %5.3f , RMS : %5.3f \n", fhPt->GetMean(), fhPt->GetRMS() ) ;
6639984f	326
	327	TCanvas * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
	328	c->Divide(1, 3);
	329
	330	c->cd(1) ;
	331	gPad->SetLogy();
	332	fhPt->SetLineColor(2);
	333	fhPt->Draw();
	334
	335	c->cd(2) ;
	336	fhPhi->SetLineColor(2);
	337	fhPhi->Draw();
	338
	339	c->cd(3) ;
	340	fhEta->SetLineColor(2);
	341	fhEta->Draw();
	342
	343	c->Print("Example.eps");
f8006433	344
	345	const Int_t buffersize = 1024;
	346	char line[buffersize] ;
	347	snprintf(line,buffersize, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
6639984f	348	gROOT->ProcessLine(line);
f8006433	349	snprintf(line,buffersize, ".!rm -fR *.eps");
6639984f	350	gROOT->ProcessLine(line);
6639984f	351
a3aebfff	352	printf("AliAnaExample::Terminate() - !! All the eps files are in %s.tar.gz !!!", GetName());
6639984f	353
6639984f	354	}