1 - Rename AliAnaMaker and AliAnaBaseClass to AliAnaPartCorrMaker and AliAnaPartCorrB...
[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 //_________________________________________________________________________
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 //-- Author: Gustavo Conesa (LNF-INFN) 
27 //_________________________________________________________________________
28
29
30 // --- ROOT system ---
31 #include "Riostream.h"
32 #include "TClonesArray.h"
33 #include "TH2F.h"
34 #include "TParticle.h"
35
36 //---- AliRoot system ----
37 #include "AliAnaExample.h"
38 #include "AliCaloTrackReader.h"
39 #include "AliLog.h"
40 #include "AliAODParticleCorrelation.h"
41 #include "AliESDCaloCells.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   fhPt  = new TH1F ("hPt","p_T distribution", 100,0,100); 
109   fhPt->SetXTitle("p_{T} (GeV/c)");
110   outputContainer->Add(fhPt);
111  
112   fhPhi  = new TH1F ("hPhi","#phi distribution", 100,0,TMath::TwoPi()); 
113   fhPhi->SetXTitle("#phi (rad)");
114   outputContainer->Add(fhPhi);
115
116   fhEta  = new TH1F ("hEta","#eta distribution", 100,-2,2); 
117   fhEta->SetXTitle("#eta ");
118   outputContainer->Add(fhEta);
119
120   if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
121     //Calo cells
122     fhNCells  = new TH1F ("hNCells","# cells per event", 100,0,1000); 
123     fhNCells->SetXTitle("n cells");
124     outputContainer->Add(fhNCells);
125     
126     fhAmplitude  = new TH1F ("hAmplitude","#eta distribution", 100,0,1000); 
127     fhAmplitude->SetXTitle("Amplitude ");
128     outputContainer->Add(fhAmplitude);
129   } 
130   
131   if(IsDataMC()){
132     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", 100,0,100,100,0,100); 
133     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
134     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
135     outputContainer->Add(fh2Pt);
136     
137     fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", 100,0,TMath::TwoPi(), 100,0,TMath::TwoPi()); 
138     fh2Phi->SetXTitle("#phi_{rec} (rad)");
139     fh2Phi->SetYTitle("#phi_{gen} (rad)");
140     outputContainer->Add(fh2Phi);
141     
142     fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", 100,-1,1,100,-1,1); 
143     fh2Eta->SetXTitle("#eta_{rec} ");
144     fh2Eta->SetYTitle("#eta_{gen} ");
145     outputContainer->Add(fh2Eta);
146     
147   }
148   return outputContainer;
149 }
150
151  //__________________________________________________
152 void AliAnaExample::InitParameters()
153
154   //Initialize the parameters of the analysis.
155   
156   fPdg = 22; //Keep photons
157   fDetector = "PHOS";
158   
159 }
160
161 //__________________________________________________________________
162 void AliAnaExample::Print(const Option_t * opt) const
163 {
164   //Print some relevant parameters set for the analysis
165   if(! opt)
166     return;
167   
168   printf("Min Pt = %3.2f\n", GetMinPt());
169   printf("Max Pt = %3.2f\n", GetMaxPt());
170   printf("Select clusters with pdg %d \n",fPdg);
171   printf("Select detector %s \n",fDetector.Data());
172   
173
174
175 //__________________________________________________________________
176 void  AliAnaExample::MakeAnalysisFillAOD() 
177 {
178   //Do analysis and fill aods
179   
180   //Some prints
181   if(GetDebug() > 0){
182     if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in emcal aod entries %d\n", GetAODEMCAL()->GetEntries());
183     if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntries());
184     if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntries());
185   }
186   
187   //Get List with tracks or clusters  
188   TClonesArray * partList = new TClonesArray;
189   if(fDetector == "CTS") partList = GetAODCTS();
190   else if(fDetector == "EMCAL") partList = GetAODEMCAL();
191   else if(fDetector == "PHOS") partList = GetAODPHOS();
192   
193   if(!partList || partList->GetEntries() == 0) return ;
194   
195   //Fill AODCaloClusters and AODParticleCorrelation with PHOS/EMCAL aods
196   if(fDetector == "EMCAL" || fDetector == "PHOS"){
197     
198     //WORK WITH CALOCLUSTERS
199     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
200     ConnectAODCaloClusters(); //Do Only when filling AODCaloClusters 
201     if(GetDebug() > 0) printf("Example: in calo clusters aod entries %d\n", GetAODCaloClusters()->GetEntries());
202     
203     //Get vertex for photon momentum calculation
204     Double_t v[3] ; //vertex ;
205     GetReader()->GetVertex(v);
206     
207     TLorentzVector mom ;
208     for(Int_t i = 0; i < partList->GetEntries(); i++){
209       
210       AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (partList->At(i));
211       
212       //Fill AODCaloClusters  
213       if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
214         AddAODCaloCluster(AliAODCaloCluster(*(calo)));
215       
216       //Fill AODParticleCorrelation after some selection
217       calo->GetMomentum(mom,v);
218       Int_t pdg = fPdg;
219
220       if(IsCaloPIDOn()){
221         Double_t pid[13];
222         calo->GetPID(pid);
223         pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
224         cout<<"PDG "<<pdg<<endl;
225         //pdg = GetCaloPID()->GetPdg(fDetector,mom,
226         //                calo->GetM02(), calo->GetM02(),
227         //                calo->GetDispersion(), 0, 0); 
228       }
229       
230       //Acceptance selection   
231       Bool_t in = kTRUE;
232       if(IsFidutialCutOn())
233         in =  GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
234
235       if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
236
237       //Select cluster if momentum, pdg and acceptance are good
238       if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
239         AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
240         //AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
241         ph.SetLabel(calo->GetLabel(0));
242         ph.SetPdg(pdg);
243         ph.SetDetector(fDetector);
244         AddAODParticleCorrelation(ph);
245       }//selection
246     }//loop
247     
248     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
249     //WORK WITH ESDCALOCELLS
250     //Don't connect in the same analysis PHOS and EMCAL cells.
251
252     AliESDCaloCells * esdCell = new AliESDCaloCells ;
253     if(fDetector == "PHOS") {
254       ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
255       esdCell = (AliESDCaloCells *) GetPHOSCells();
256     }
257     else  {
258       ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
259       esdCell = (AliESDCaloCells *) GetEMCALCells();
260     }
261     //Some prints
262     if(GetDebug() > 0 && esdCell )
263       printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());    
264     
265     //Fill AODCells in file
266     Int_t ncells = esdCell->GetNumberOfCells() ;
267     GetAODCaloCells()->CreateContainer(ncells);
268
269     GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
270
271     for (Int_t iCell = 0; iCell < ncells; iCell++) {      
272       if(GetDebug() > 2)  printf("cell : amp %f, absId %d,  time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
273
274       GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
275     }
276     GetAODCaloCells()->Sort();
277     } 
278   }//cluster-cell analysis
279   else if(fDetector == "CTS"){ //Track analysis
280     //Fill AODParticleCorrelation with CTS aods
281     TVector3 p3;
282     for(Int_t i = 0; i < GetAODCTS()->GetEntries(); i++){
283       
284       AliAODTrack * track =  dynamic_cast<AliAODTrack*> (GetAODCTS()->At(i));
285       
286       //Fill AODParticleCorrelation after some selection       
287       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
288       p3.SetXYZ(mom[0],mom[1],mom[2]);
289       
290       //Acceptance selection
291       Bool_t in =  GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
292       if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
293       if(p3.Pt() > GetMinPt() && in) {
294         AliAODParticleCorrelation tr = AliAODParticleCorrelation(mom[0],mom[1],mom[2],0);
295         //AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
296         tr.SetDetector("CTS");
297         AddAODParticleCorrelation(tr);
298         
299       }//selection
300     }//loop
301   }//CTS analysis
302   
303   if(GetDebug() > 0) {
304     printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntries());
305     printf("Example: final aod branch entries %d\n", GetAODBranch()->GetEntries());
306     printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
307   }
308
309
310 //__________________________________________________________________
311 void  AliAnaExample::MakeAnalysisFillHistograms() 
312 {
313   //Do analysis and fill histograms
314
315   //Loop on stored AODParticles
316   Int_t naod = GetAODBranch()->GetEntriesFast();
317   if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
318   for(Int_t iaod = 0; iaod < naod ; iaod++){
319    AliAODParticleCorrelation* ph =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
320
321    fhPt->Fill(ph->Pt());
322    fhPhi->Fill(ph->Phi());
323    fhEta->Fill(ph->Eta());
324
325    if(IsDataMC()){
326      //Play with the MC stack if available
327      AliStack * stack =  GetMCStack() ;
328      
329      if(ph->GetLabel() < 0 || !stack) {
330        printf("*** bad label or no stack ***:  label %d \n", ph->GetLabel());
331        continue;
332      }
333      
334      if(ph->GetLabel() >=  stack->GetNtrack()) {
335        printf("*** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
336        continue ;
337      }
338      
339      TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
340      
341      fh2Pt->Fill(ph->Pt(), mom->Pt());
342      fh2Phi->Fill(ph->Phi(), mom->Phi());
343      fh2Eta->Fill(ph->Eta(), mom->Eta());
344    }//Work with stack also
345   }// aod branch loop
346
347   // CaloCells histograms
348   if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
349     if(GetAODCaloCells()){
350       
351       Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
352       fhNCells->Fill(ncells) ;
353       
354       for (Int_t iCell = 0; iCell < ncells; iCell++) {      
355         if(GetDebug() > 2)  printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
356         fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
357       }
358     }//calo cells container exist
359   }
360 }