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