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