Correct and clean the vertex retrieval in case of SE or ME analysis
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaExample.cxx
CommitLineData
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"
43#include "AliCaloPID.h"
ff45398a 44#include "AliFiducialCut.h"
f37fa8d2 45#include "AliVCluster.h"
477d6cee 46#include "AliAODTrack.h"
1c5acb87 47
48ClassImp(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//________________________________________________________________________
61TList * 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}
111
477d6cee 112//__________________________________________________
1c5acb87 113void 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//__________________________________________________________________
126void 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//__________________________________________________________________
141void AliAnaExample::MakeAnalysisFillAOD()
142{
477d6cee 143 //Do analysis and fill aods
144
145 //Some prints
146 if(GetDebug() > 0){
a3aebfff 147 if(fDetector == "EMCAL" && GetAODEMCAL())printf("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
f37fa8d2 148 if(fDetector == "CTS" && GetAODCTS()) printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS() ->GetEntriesFast());
149 if(fDetector == "PHOS" && GetAODPHOS()) printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS() ->GetEntriesFast());
477d6cee 150 }
151
152 //Get List with tracks or clusters
4a745797 153 TObjArray * partList = 0x0;
f37fa8d2 154 if(fDetector == "CTS") partList = GetAODCTS();
477d6cee 155 else if(fDetector == "EMCAL") partList = GetAODEMCAL();
f37fa8d2 156 else if(fDetector == "PHOS") partList = GetAODPHOS();
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
211 }//loop
f37fa8d2 212 }//Calorimeters
477d6cee 213 else if(fDetector == "CTS"){ //Track analysis
214 //Fill AODParticle with CTS aods
215 TVector3 p3;
216 for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
217
218 AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(i));
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") ;
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//__________________________________________________________________
244void 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//________________________________________________________________________
293void 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 316void 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()) ;
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);
351
a3aebfff 352 printf("AliAnaExample::Terminate() - !! All the eps files are in %s.tar.gz !!!", GetName());
6639984f 353
354}