[u/mrichter/AliRoot.git] / PWG4 / PartCorr / 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)
//-- Author: Gustavo Conesa (LNF-INFN) 
//_________________________________________________________________________


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

//---- AliRoot system ----
#include "AliAnaExample.h"
#include "AliCaloTrackReader.h"
#include "AliLog.h"
#include "AliAODParticleCorrelation.h"
#include "AliESDCaloCells.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") ; 

  fhPt  = new TH1F ("hPt","p_T distribution", 100,0,100); 
  fhPt->SetXTitle("p_{T} (GeV/c)");
  outputContainer->Add(fhPt);
 
  fhPhi  = new TH1F ("hPhi","#phi distribution", 100,0,TMath::TwoPi()); 
  fhPhi->SetXTitle("#phi (rad)");
  outputContainer->Add(fhPhi);

  fhEta  = new TH1F ("hEta","#eta distribution", 100,-2,2); 
  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", 100,0,100,100,0,100); 
    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", 100,0,TMath::TwoPi(), 100,0,TMath::TwoPi()); 
    fh2Phi->SetXTitle("#phi_{rec} (rad)");
    fh2Phi->SetYTitle("#phi_{gen} (rad)");
    outputContainer->Add(fh2Phi);
    
    fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", 100,-1,1,100,-1,1); 
    fh2Eta->SetXTitle("#eta_{rec} ");
    fh2Eta->SetYTitle("#eta_{gen} ");
    outputContainer->Add(fh2Eta);
    
  }
  return outputContainer;
}

 //__________________________________________________
void AliAnaExample::InitParameters()
{ 
  //Initialize the parameters of the analysis.
  
  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 AODParticleCorrelation 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 =  dynamic_cast<AliAODCaloCluster*> (partList->At(i));
      
      //Fill AODCaloClusters  
      if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
	AddAODCaloCluster(AliAODCaloCluster(*(calo)));
      
      //Fill AODParticleCorrelation 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());
	cout<<"PDG "<<pdg<<endl;
	//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) {
	AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
	//AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
	ph.SetLabel(calo->GetLabel(0));
	ph.SetPdg(pdg);
	ph.SetDetector(fDetector);
	AddAODParticleCorrelation(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();
    }
    //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 AODParticleCorrelation with CTS aods
    TVector3 p3;
    for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
      
      AliAODTrack * track =  dynamic_cast<AliAODTrack*> (GetAODCTS()->At(i));
      
      //Fill AODParticleCorrelation 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) {
	AliAODParticleCorrelation tr = AliAODParticleCorrelation(mom[0],mom[1],mom[2],0);
	//AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
	tr.SetDetector("CTS");
	AddAODParticleCorrelation(tr);
	
      }//selection
    }//loop
  }//CTS analysis
  
  if(GetDebug() > 0) {
    printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
    printf("Example: final aod branch entries %d\n", GetAODBranch()->GetEntriesFast());
    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 = GetAODBranch()->GetEntriesFast();
  if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
  for(Int_t iaod = 0; iaod < naod ; iaod++){
   AliAODParticleCorrelation* ph =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->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
f3050824	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	//-- Author: Gustavo Conesa (LNF-INFN)
	27	//_________________________________________________________________________
	28
	29
	30	// --- ROOT system ---
	31	#include "Riostream.h"
	32	#include "TClonesArray.h"
	33	#include "TH2F.h"
	34	#include "TParticle.h"
	35
	36	//---- AliRoot system ----
	37	#include "AliAnaExample.h"
	38	#include "AliCaloTrackReader.h"
	39	#include "AliLog.h"
	40	#include "AliAODParticleCorrelation.h"
	41	#include "AliESDCaloCells.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	fhPt = new TH1F ("hPt","p_T distribution", 100,0,100);
109	fhPt->SetXTitle("p_{T} (GeV/c)");
110	outputContainer->Add(fhPt);
111
112	fhPhi = new TH1F ("hPhi","#phi distribution", 100,0,TMath::TwoPi());
113	fhPhi->SetXTitle("#phi (rad)");
114	outputContainer->Add(fhPhi);
115
116	fhEta = new TH1F ("hEta","#eta distribution", 100,-2,2);
117	fhEta->SetXTitle("#eta ");
118	outputContainer->Add(fhEta);
119
120	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
121	//Calo cells
122	fhNCells = new TH1F ("hNCells","# cells per event", 100,0,1000);
123	fhNCells->SetXTitle("n cells");
124	outputContainer->Add(fhNCells);
125
126	fhAmplitude = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
127	fhAmplitude->SetXTitle("Amplitude ");
128	outputContainer->Add(fhAmplitude);
129	}
130
131	if(IsDataMC()){
132	fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", 100,0,100,100,0,100);
133	fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
134	fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
135	outputContainer->Add(fh2Pt);
136
137	fh2Phi = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", 100,0,TMath::TwoPi(), 100,0,TMath::TwoPi());
138	fh2Phi->SetXTitle("#phi_{rec} (rad)");
139	fh2Phi->SetYTitle("#phi_{gen} (rad)");
140	outputContainer->Add(fh2Phi);
141
142	fh2Eta = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", 100,-1,1,100,-1,1);
143	fh2Eta->SetXTitle("#eta_{rec} ");
144	fh2Eta->SetYTitle("#eta_{gen} ");
145	outputContainer->Add(fh2Eta);
146
147	}
148	return outputContainer;
149	}
150
151	//__________________________________________________
152	void AliAnaExample::InitParameters()
153	{
154	//Initialize the parameters of the analysis.
155
156	fPdg = 22; //Keep photons
157	fDetector = "PHOS";
158
159	}
160
161	//__________________________________________________________________
162	void AliAnaExample::Print(const Option_t * opt) const
163	{
164	//Print some relevant parameters set for the analysis
165	if(! opt)
166	return;
167
168	printf("Min Pt = %3.2f\n", GetMinPt());
169	printf("Max Pt = %3.2f\n", GetMaxPt());
170	printf("Select clusters with pdg %d \n",fPdg);
171	printf("Select detector %s \n",fDetector.Data());
172
173	}
174
175	//__________________________________________________________________
176	void AliAnaExample::MakeAnalysisFillAOD()
177	{
178	//Do analysis and fill aods
179
180	//Some prints
181	if(GetDebug() > 0){
182	if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in emcal aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
183	if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
184	if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
185	}
186
187	//Get List with tracks or clusters
188	TClonesArray * partList = new TClonesArray;
189	if(fDetector == "CTS") partList = GetAODCTS();
190	else if(fDetector == "EMCAL") partList = GetAODEMCAL();
191	else if(fDetector == "PHOS") partList = GetAODPHOS();
192
193	if(!partList \|\| partList->GetEntriesFast() == 0) return ;
194
195	//Fill AODCaloClusters and AODParticleCorrelation with PHOS/EMCAL aods
196	if(fDetector == "EMCAL" \|\| fDetector == "PHOS"){
197
198	//WORK WITH CALOCLUSTERS
199	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
200	ConnectAODCaloClusters(); //Do Only when filling AODCaloClusters
201	if(GetDebug() > 0) printf("Example: in calo clusters aod entries %d\n", GetAODCaloClusters()->GetEntriesFast());
202
203	//Get vertex for photon momentum calculation
204	Double_t v[3] ; //vertex ;
205	GetReader()->GetVertex(v);
206
207	TLorentzVector mom ;
208	for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
209
210	AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (partList->At(i));
211
212	//Fill AODCaloClusters
213	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
214	AddAODCaloCluster(AliAODCaloCluster(*(calo)));
215
216	//Fill AODParticleCorrelation after some selection
217	calo->GetMomentum(mom,v);
218	Int_t pdg = fPdg;
219
220	if(IsCaloPIDOn()){
221	Double_t pid[13];
222	calo->GetPID(pid);
223	pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
224	cout<<"PDG "<<pdg<<endl;
225	//pdg = GetCaloPID()->GetPdg(fDetector,mom,
226	// calo->GetM02(), calo->GetM02(),
227	// calo->GetDispersion(), 0, 0);
228	}
229
230	//Acceptance selection
231	Bool_t in = kTRUE;
232	if(IsFidutialCutOn())
233	in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
234
235	if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
236
237	//Select cluster if momentum, pdg and acceptance are good
238	if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
239	AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
240	//AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
241	ph.SetLabel(calo->GetLabel(0));
242	ph.SetPdg(pdg);
243	ph.SetDetector(fDetector);
244	AddAODParticleCorrelation(ph);
245	}//selection
246	}//loop
247
248	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
249	//WORK WITH ESDCALOCELLS
250	//Don't connect in the same analysis PHOS and EMCAL cells.
251
252	AliESDCaloCells * esdCell = new AliESDCaloCells ;
253	if(fDetector == "PHOS") {
254	ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
255	esdCell = (AliESDCaloCells *) GetPHOSCells();
256	}
257	else {
258	ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
259	esdCell = (AliESDCaloCells *) GetEMCALCells();
260	}
261	//Some prints
262	if(GetDebug() > 0 && esdCell )
263	printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
264
265	//Fill AODCells in file
266	Int_t ncells = esdCell->GetNumberOfCells() ;
267	GetAODCaloCells()->CreateContainer(ncells);
268
269	GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
270
271	for (Int_t iCell = 0; iCell < ncells; iCell++) {
272	if(GetDebug() > 2) printf("cell : amp %f, absId %d, time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
273
274	GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
275	}
276	GetAODCaloCells()->Sort();
277	}
278	}//cluster-cell analysis
279	else if(fDetector == "CTS"){ //Track analysis
280	//Fill AODParticleCorrelation with CTS aods
281	TVector3 p3;
282	for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
283
284	AliAODTrack * track = dynamic_cast<AliAODTrack*> (GetAODCTS()->At(i));
285
286	//Fill AODParticleCorrelation after some selection
287	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
288	p3.SetXYZ(mom[0],mom[1],mom[2]);
289
290	//Acceptance selection
291	Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
292	if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
293	if(p3.Pt() > GetMinPt() && in) {
294	AliAODParticleCorrelation tr = AliAODParticleCorrelation(mom[0],mom[1],mom[2],0);
295	//AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
296	tr.SetDetector("CTS");
297	AddAODParticleCorrelation(tr);
298
299	}//selection
300	}//loop
301	}//CTS analysis
302
303	if(GetDebug() > 0) {
304	printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
305	printf("Example: final aod branch entries %d\n", GetAODBranch()->GetEntriesFast());
306	printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
307	}
308	}
309
310	//__________________________________________________________________
311	void AliAnaExample::MakeAnalysisFillHistograms()
312	{
313	//Do analysis and fill histograms
314
315	//Loop on stored AODParticles
316	Int_t naod = GetAODBranch()->GetEntriesFast();
317	if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
318	for(Int_t iaod = 0; iaod < naod ; iaod++){
319	AliAODParticleCorrelation* ph = dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
320
321	fhPt->Fill(ph->Pt());
322	fhPhi->Fill(ph->Phi());
323	fhEta->Fill(ph->Eta());
324
325	if(IsDataMC()){
326	//Play with the MC stack if available
327	AliStack * stack = GetMCStack() ;
328
329	if(ph->GetLabel() < 0 \|\| !stack) {
330	printf("* bad label or no stack *: label %d \n", ph->GetLabel());
331	continue;
332	}
333
334	if(ph->GetLabel() >= stack->GetNtrack()) {
335	printf("* large label *: label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
336	continue ;
337	}
338
339	TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
340
341	fh2Pt->Fill(ph->Pt(), mom->Pt());
342	fh2Phi->Fill(ph->Phi(), mom->Phi());
343	fh2Eta->Fill(ph->Eta(), mom->Eta());
344	}//Work with stack also
345	}// aod branch loop
346
347	// CaloCells histograms
348	if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
349	if(GetAODCaloCells()){
350
351	Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
352	fhNCells->Fill(ncells) ;
353
354	for (Int_t iCell = 0; iCell < ncells; iCell++) {
355	if(GetDebug() > 2) printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
356	fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
357	}
358	}//calo cells container exist
359	}
360	}