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 "TRefArray.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 "AliFidutialCut.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 = GetHistoNPtBins();
114 Int_t nphibins = GetHistoNPhiBins();
115 Int_t netabins = GetHistoNEtaBins();
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);
164 return outputContainer;
167 //__________________________________________________
168 void AliAnaExample::InitParameters()
170 //Initialize the parameters of the analysis.
171 SetOutputAODClassName("AliAODPWG4Particle");
172 SetOutputAODName("example");
173 fPdg = 22; //Keep photons
178 //__________________________________________________________________
179 void AliAnaExample::Print(const Option_t * opt) const
181 //Print some relevant parameters set for the analysis
185 printf("Min Pt = %3.2f\n", GetMinPt());
186 printf("Max Pt = %3.2f\n", GetMaxPt());
187 printf("Select clusters with pdg %d \n",fPdg);
188 printf("Select detector %s \n",fDetector.Data());
192 //__________________________________________________________________
193 void AliAnaExample::MakeAnalysisFillAOD()
195 //Do analysis and fill aods
199 if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
200 if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
201 if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
204 //Get List with tracks or clusters
205 TRefArray * partList = new TRefArray();
206 if(fDetector == "CTS") partList = GetAODCTS();
207 else if(fDetector == "EMCAL") partList = GetAODEMCAL();
208 else if(fDetector == "PHOS") partList = GetAODPHOS();
210 if(!partList || partList->GetEntriesFast() == 0) return ;
212 //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
213 if(fDetector == "EMCAL" || fDetector == "PHOS"){
215 //Get vertex for photon momentum calculation
216 Double_t v[3] ; //vertex ;
217 GetReader()->GetVertex(v);
220 for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
222 AliAODCaloCluster * calo = (AliAODCaloCluster*) (partList->At(i));
224 //Fill AODParticle after some selection
225 calo->GetMomentum(mom,v);
231 pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
232 //pdg = GetCaloPID()->GetPdg(fDetector,mom,
233 // calo->GetM02(), calo->GetM02(),
234 // calo->GetDispersion(), 0, 0);
237 //Acceptance selection
239 if(IsFidutialCutOn())
240 in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
242 if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
244 //Select cluster if momentum, pdg and acceptance are good
245 if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
246 AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
247 //AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
248 ph.SetLabel(calo->GetLabel(0));
250 ph.SetDetector(fDetector);
255 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
256 //WORK WITH ESDCALOCELLS
257 //Don't connect in the same analysis PHOS and EMCAL cells.
259 AliESDCaloCells * esdCell = new AliESDCaloCells ;
260 if(fDetector == "PHOS") {
261 //ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
262 esdCell = (AliESDCaloCells *) GetPHOSCells();
265 //ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
266 esdCell = (AliESDCaloCells *) GetEMCALCells();
270 printf("FillAOD ABORT: No CELLS available for analysis");
274 if(GetDebug() > 0 && esdCell )
275 printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
277 //Fill AODCells in file
278 Int_t ncells = esdCell->GetNumberOfCells() ;
279 GetAODCaloCells()->CreateContainer(ncells);
281 GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
283 for (Int_t iCell = 0; iCell < ncells; iCell++) {
284 if(GetDebug() > 2) printf("cell : amp %f, absId %d, time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
286 GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
288 GetAODCaloCells()->Sort();
290 }//cluster-cell analysis
291 else if(fDetector == "CTS"){ //Track analysis
292 //Fill AODParticle with CTS aods
294 for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
296 AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(i));
298 //Fill AODParticle after some selection
299 Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
300 p3.SetXYZ(mom[0],mom[1],mom[2]);
302 //Acceptance selection
303 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
304 if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
305 if(p3.Pt() > GetMinPt() && in) {
306 AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
307 tr.SetDetector("CTS");
314 if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
315 //printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
316 printf("Example: final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
317 // if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
318 //printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
322 //__________________________________________________________________
323 void AliAnaExample::MakeAnalysisFillHistograms()
325 //Do analysis and fill histograms
327 //Loop on stored AODParticles
328 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
329 if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
330 for(Int_t iaod = 0; iaod < naod ; iaod++){
331 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
333 fhPt->Fill(ph->Pt());
334 fhPhi->Fill(ph->Phi());
335 fhEta->Fill(ph->Eta());
338 //Play with the MC stack if available
339 AliStack * stack = GetMCStack() ;
341 if(ph->GetLabel() < 0 || !stack) {
342 printf("*** bad label or no stack ***: label %d \n", ph->GetLabel());
346 if(ph->GetLabel() >= stack->GetNtrack()) {
347 printf("*** large label ***: label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
351 TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
353 fh2Pt->Fill(ph->Pt(), mom->Pt());
354 fh2Phi->Fill(ph->Phi(), mom->Phi());
355 fh2Eta->Fill(ph->Eta(), mom->Eta());
356 }//Work with stack also
359 // CaloCells histograms
360 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
361 if(GetAODCaloCells()){
363 Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
364 fhNCells->Fill(ncells) ;
366 for (Int_t iCell = 0; iCell < ncells; iCell++) {
367 if(GetDebug() > 2) printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
368 fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
370 }//calo cells container exist
374 //__________________________________________________________________
375 void AliAnaExample::Terminate()
378 //Do some plots to end
380 printf(" *** %s Report:", GetName()) ;
381 printf(" pt : %5.3f , RMS : %5.3f \n", fhPt->GetMean(), fhPt->GetRMS() ) ;
383 TCanvas * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
388 fhPt->SetLineColor(2);
392 fhPhi->SetLineColor(2);
396 fhEta->SetLineColor(2);
399 c->Print("Example.eps");
402 sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
403 gROOT->ProcessLine(line);
404 sprintf(line, ".!rm -fR *.eps");
405 gROOT->ProcessLine(line);
407 printf("!! All the eps files are in %s.tar.gz !!!", GetName());