1)Terminate() method implemented in the frame. Simple examples on what to do with...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaExample.cxx
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"
33 #include "TCanvas.h"
34 #include "TPad.h"
35 #include "TROOT.h"
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
46 ClassImp(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 //____________________________________________________________________________
60 AliAnaExample::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 //_________________________________________________________________________
69 AliAnaExample & 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 //________________________________________________________________________
101 TList *  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  //__________________________________________________
166 void 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 //__________________________________________________________________
177 void 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 //__________________________________________________________________
191 void  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());  
322                 //      if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
323                 //printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
324         }
325
326
327 //__________________________________________________________________
328 void  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 }
378
379 //__________________________________________________________________
380 void  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 }