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