1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
25 //-- Author: Gustavo Conesa (INFN-LNF)
26 //_________________________________________________________________________
29 // --- ROOT system ---
30 //#include "Riostream.h"
31 #include "TObjArray.h"
32 #include "TParticle.h"
38 //---- AliRoot system ----
39 #include "AliAnaExample.h"
40 #include "AliCaloTrackReader.h"
41 #include "AliAODPWG4Particle.h"
42 #include "AliESDCaloCells.h"
44 #include "AliCaloPID.h"
45 #include "AliFiducialCut.h"
46 #include "AliAODCaloCells.h"
47 #include "AliAODCaloCluster.h"
48 #include "AliAODTrack.h"
50 ClassImp(AliAnaExample)
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)
59 //Initialize parameters
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)
72 //_________________________________________________________________________
73 AliAnaExample & AliAnaExample::operator = (const AliAnaExample & ex)
75 // assignment operator
77 if(this == &ex)return *this;
78 ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
81 fDetector = ex.fDetector;
88 fhNCells = ex.fhNCells;
89 fhAmplitude = ex.fhAmplitude;
95 // //____________________________________________________________________________
96 // AliAnaExample::~AliAnaExample()
98 // // Remove all pointers except analysis output pointers.
104 //________________________________________________________________________
105 TList * AliAnaExample::GetCreateOutputObjects()
107 // Create histograms to be saved in output file and
108 // store them in fOutputContainer
110 TList * outputContainer = new TList() ;
111 outputContainer->SetName("ExampleHistos") ;
113 Int_t nptbins = GetHistoPtBins();
114 Int_t nphibins = GetHistoPhiBins();
115 Int_t netabins = GetHistoEtaBins();
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();
123 fhPt = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
124 fhPt->SetXTitle("p_{T} (GeV/c)");
125 outputContainer->Add(fhPt);
127 fhPhi = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
128 fhPhi->SetXTitle("#phi (rad)");
129 outputContainer->Add(fhPhi);
131 fhEta = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
132 fhEta->SetXTitle("#eta ");
133 outputContainer->Add(fhEta);
135 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
138 fhNCells = new TH1F ("hNCells","# cells per event", 100,0,1000);
139 fhNCells->SetXTitle("n cells");
140 outputContainer->Add(fhNCells);
142 fhAmplitude = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
143 fhAmplitude->SetXTitle("Amplitude ");
144 outputContainer->Add(fhAmplitude);
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);
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);
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);
165 return outputContainer;
168 //__________________________________________________
169 void AliAnaExample::InitParameters()
171 //Initialize the parameters of the analysis.
172 SetOutputAODClassName("AliAODPWG4Particle");
173 SetOutputAODName("PWG4Particle");
174 AddToHistogramsName("AnaExample_");
176 fPdg = 22; //Keep photons
181 //__________________________________________________________________
182 void AliAnaExample::Print(const Option_t * opt) const
184 //Print some relevant parameters set for the analysis
188 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
189 AliAnaPartCorrBaseClass::Print(" ");
191 printf("Select clusters with pdg %d \n",fPdg);
192 printf("Select detector %s \n",fDetector.Data());
196 //__________________________________________________________________
197 void AliAnaExample::MakeAnalysisFillAOD()
199 //Do analysis and fill aods
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());
208 //Get List with tracks or clusters
209 TObjArray * partList = 0x0;
210 if(fDetector == "CTS") partList = GetAODCTS();
211 else if(fDetector == "EMCAL") partList = GetAODEMCAL();
212 else if(fDetector == "PHOS") partList = GetAODPHOS();
214 if(!partList || partList->GetEntriesFast() == 0) return ;
216 //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
217 if(fDetector == "EMCAL" || fDetector == "PHOS"){
219 //Get vertex for photon momentum calculation
220 Double_t v[3] ; //vertex ;
221 GetReader()->GetVertex(v);
224 for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
226 AliAODCaloCluster * calo = (AliAODCaloCluster*) (partList->At(i));
228 //Fill AODParticle after some selection
229 calo->GetMomentum(mom,v);
233 const Double_t *pid = calo->GetPID();
234 pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
235 //pdg = GetCaloPID()->GetPdg(fDetector,mom,
236 // calo->GetM02(), calo->GetM02(),
237 // calo->GetDispersion(), 0, 0);
240 //Acceptance selection
242 if(IsFiducialCutOn())
243 in = GetFiducialCut()->IsInFiducialCut(mom,fDetector) ;
245 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);
247 //Select cluster if momentum, pdg and acceptance are good
248 if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
249 AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
250 //AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
251 ph.SetLabel(calo->GetLabel());
253 ph.SetDetector(fDetector);
258 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
259 //WORK WITH ESDCALOCELLS
260 //Don't connect in the same analysis PHOS and EMCAL cells.
262 AliESDCaloCells * esdCell = new AliESDCaloCells ;
263 if(fDetector == "PHOS") {
264 ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
265 esdCell = (AliESDCaloCells *) GetPHOSCells();
268 ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
269 esdCell = (AliESDCaloCells *) GetEMCALCells();
273 printf("AliAnaExample::MakeAnalysisFillAOD() - STOP: No CELLS available for analysis");
277 if(GetDebug() > 0 && esdCell )
278 printf("AliAnaExample::MakeAnalysisFillAOD() - In ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
280 //Fill AODCells in file
281 Int_t ncells = esdCell->GetNumberOfCells() ;
282 GetAODCaloCells()->CreateContainer(ncells);
284 GetAODCaloCells()->SetType(esdCell->GetType());
286 for (Int_t iCell = 0; iCell < ncells; iCell++) {
287 if(GetDebug() > 2) printf("AliAnaExample::MakeAnalysisFillAOD() - Cell : amp %f, absId %d, time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
289 GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
291 GetAODCaloCells()->Sort();
293 }//cluster-cell analysis
294 else if(fDetector == "CTS"){ //Track analysis
295 //Fill AODParticle with CTS aods
297 for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
299 AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(i));
301 //Fill AODParticle after some selection
302 Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
303 p3.SetXYZ(mom[0],mom[1],mom[2]);
305 //Acceptance selection
306 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
307 if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
308 if(p3.Pt() > GetMinPt() && in) {
309 AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
310 tr.SetDetector("CTS");
317 if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
318 //printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
319 printf("AliAnaExample::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
320 // if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
321 //printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
325 //__________________________________________________________________
326 void AliAnaExample::MakeAnalysisFillHistograms()
328 //Do analysis and fill histograms
330 //Loop on stored AODParticles
331 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
332 if(GetDebug() > 0) printf("AliAnaExample::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
333 for(Int_t iaod = 0; iaod < naod ; iaod++){
334 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
336 fhPt->Fill(ph->Pt());
337 fhPhi->Fill(ph->Phi());
338 fhEta->Fill(ph->Eta());
341 //Play with the MC stack if available
345 if(GetReader()->ReadStack()){
346 AliStack * stack = GetMCStack() ;
348 if(ph->GetLabel() < 0 || !stack) {
349 printf("AliAnaExample::MakeAnalysisFillHistograms() *** bad label or no stack ***: label %d \n", ph->GetLabel());
353 if(ph->GetLabel() >= stack->GetNtrack()) {
354 printf("AliAnaExample::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
358 TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
360 phiprim = mom->Phi();
361 etaprim = mom->Eta();
363 else if(GetReader()->ReadAODMCParticles()){
364 printf("AliAnaExample::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
366 fh2Pt->Fill(ph->Pt(), ptprim);
367 fh2Phi->Fill(ph->Phi(), phiprim);
368 fh2Eta->Fill(ph->Eta(), etaprim);
369 }//Work with stack also
372 // CaloCells histograms
373 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
374 if(GetAODCaloCells()){
376 Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
377 fhNCells->Fill(ncells) ;
379 for (Int_t iCell = 0; iCell < ncells; iCell++) {
380 if(GetDebug() > 2) printf("AliAnaExample::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
381 fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
383 }//calo cells container exist
387 //________________________________________________________________________
388 void AliAnaExample::ReadHistograms(TList* outputList)
390 // Needed when Terminate is executed in distributed environment
391 // Refill analysis histograms of this class with corresponding histograms in output list.
393 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
394 // the first one and then add the next
395 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hPt"));
397 //Read histograms, must be in the same order as in GetCreateOutputObject.
398 fhPt = (TH1F *) outputList->At(index++);
399 fhPhi = (TH1F *) outputList->At(index++);
400 fhEta = (TH1F *) outputList->At(index++);
402 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
403 fhNCells = (TH1F *) outputList->At(index++);
404 fhAmplitude = (TH1F *) outputList->At(index++);
408 fh2Pt = (TH2F *) outputList->At(index++);
409 fh2Phi = (TH2F *) outputList->At(index++);
410 fh2Eta = (TH2F *) outputList->At(index);
415 //__________________________________________________________________
416 void AliAnaExample::Terminate(TList* outputList)
419 //Do some plots to end
421 //Recover histograms from output histograms list, needed for distributed analysis.
422 ReadHistograms(outputList);
424 printf(" AliAnaExample::Terminate() *** %s Report:", GetName()) ;
425 printf(" AliAnaExample::Terminate() pt : %5.3f , RMS : %5.3f \n", fhPt->GetMean(), fhPt->GetRMS() ) ;
427 TCanvas * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
432 fhPt->SetLineColor(2);
436 fhPhi->SetLineColor(2);
440 fhEta->SetLineColor(2);
443 c->Print("Example.eps");
446 sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
447 gROOT->ProcessLine(line);
448 sprintf(line, ".!rm -fR *.eps");
449 gROOT->ProcessLine(line);
451 printf("AliAnaExample::Terminate() - !! All the eps files are in %s.tar.gz !!!", GetName());