1)Terminate() method implemented in the frame. Simple examples on what to do with...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleHadronCorrelation.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 // Class for the analysis of particle - hadron correlations
19 // Particle (for example direct gamma) must be found in a previous analysis 
20 //-- Author: Gustavo Conesa (LNF-INFN) 
21 //////////////////////////////////////////////////////////////////////////////
22
23
24 // --- ROOT system ---
25 #include "TH2F.h"
26
27 //---- ANALYSIS system ----
28 #include "AliLog.h"
29 #include "AliNeutralMesonSelection.h" 
30 #include "AliAnaParticleHadronCorrelation.h" 
31 #include "AliCaloTrackReader.h"
32 #include "AliCaloPID.h"
33 #include "AliAODPWG4ParticleCorrelation.h"
34
35 ClassImp(AliAnaParticleHadronCorrelation)
36
37
38 //____________________________________________________________________________
39   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation() : 
40     AliAnaPartCorrBaseClass(),
41     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),  
42         fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0), 
43     fhDeltaPhiCharged(0), fhDeltaPhiNeutral(0), 
44     fhDeltaEtaCharged(0), fhDeltaEtaNeutral(0),
45     fhDeltaPhiChargedPt(0), fhDeltaPhiNeutralPt(0), 
46     fhPtImbalanceNeutral(0), fhPtImbalanceCharged(0)
47 {
48   //Default Ctor
49  
50   //Initialize parameters
51   InitParameters();
52 }
53
54 //____________________________________________________________________________
55 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(const AliAnaParticleHadronCorrelation & g) :   
56   AliAnaPartCorrBaseClass(g),
57   fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut), 
58   fSelectIsolated(g.fSelectIsolated),
59   fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral), 
60   fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral), 
61   fhDeltaPhiCharged(g.fhDeltaPhiCharged),  
62   fhDeltaPhiNeutral(g.fhDeltaPhiNeutral), 
63   fhDeltaEtaCharged(g.fhDeltaEtaCharged), 
64   fhDeltaEtaNeutral(g.fhDeltaEtaNeutral), 
65   fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt), 
66   fhDeltaPhiNeutralPt(g.fhDeltaPhiNeutralPt), 
67   fhPtImbalanceNeutral(g.fhPtImbalanceNeutral), 
68   fhPtImbalanceCharged(g.fhPtImbalanceCharged)
69 {
70   // cpy ctor
71
72 }
73
74 //_________________________________________________________________________
75 AliAnaParticleHadronCorrelation & AliAnaParticleHadronCorrelation::operator = (const AliAnaParticleHadronCorrelation & source)
76 {
77   // assignment operator
78
79   if(this == &source)return *this;
80   ((AliAnaPartCorrBaseClass *)this)->operator=(source);
81   
82   fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ; 
83   fDeltaPhiMinCut = source.fDeltaPhiMinCut ; 
84   fSelectIsolated = source.fSelectIsolated ;
85  
86   fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ; 
87   fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ; 
88   fhDeltaPhiCharged = source.fhDeltaPhiCharged ;  
89   fhDeltaPhiNeutral = source.fhDeltaPhiNeutral ; 
90   fhDeltaPhiNeutralPt = source.fhDeltaPhiNeutralPt ; 
91   fhDeltaEtaCharged = source.fhDeltaEtaCharged ; 
92   fhDeltaEtaNeutral = source.fhDeltaEtaNeutral ; 
93   fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
94
95   fhPtImbalanceNeutral = source.fhPtImbalanceNeutral ; 
96   fhPtImbalanceCharged = source.fhPtImbalanceCharged ; 
97
98   return *this;
99
100 }
101
102 //________________________________________________________________________
103 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
104 {  
105         
106         // Create histograms to be saved in output file and 
107         // store them in fOutputContainer
108         TList * outputContainer = new TList() ; 
109         outputContainer->SetName("CorrelationHistos") ; 
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         //Correlation with charged hadrons
122         if(GetReader()->IsCTSSwitchedOn()) {
123                 fhPhiCharged  = new TH2F
124                 ("PhiCharged","#phi_{h^{#pm}}  vs p_{T trigger}",
125                  nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
126                 fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
127                 fhPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
128                 
129                 fhEtaCharged  = new TH2F
130                 ("EtaCharged","#eta_{h^{#pm}}  vs p_{T trigger}",
131                  nptbins,ptmin,ptmax,netabins,etamin,etamax); 
132                 fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
133                 fhEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
134                 
135                 fhDeltaPhiCharged  = new TH2F
136                 ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
137                 nptbins,ptmin,ptmax,200,0,6.4); 
138                 fhDeltaPhiCharged->SetYTitle("#Delta #phi");
139                 fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
140                 
141                 fhDeltaPhiChargedPt  = new TH2F
142                 ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#p^{#pm}i} vs p_{T h^{#pm}}",
143                  nptbins,ptmin,ptmax,200,0,6.4);
144                 fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
145                 fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
146                 
147                 fhDeltaEtaCharged  = new TH2F
148                 ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
149                 nptbins,ptmin,ptmax,200,-2,2); 
150                 fhDeltaEtaCharged->SetYTitle("#Delta #eta");
151                 fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
152                 
153                 fhPtImbalanceCharged  = 
154                 new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
155                                 nptbins,ptmin,ptmax,1000,0.,1.2); 
156                 fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
157                 fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
158                 
159                 outputContainer->Add(fhPhiCharged) ;
160                 outputContainer->Add(fhEtaCharged) ;
161                 outputContainer->Add(fhDeltaPhiCharged) ; 
162                 outputContainer->Add(fhDeltaEtaCharged) ;
163                 outputContainer->Add(fhPtImbalanceCharged) ;
164                 outputContainer->Add(fhDeltaPhiChargedPt) ;
165         }  //Correlation with charged hadrons
166         
167         //Correlation with neutral hadrons
168         if(GetReader()->IsEMCALSwitchedOn() || GetReader()->IsPHOSSwitchedOn()){
169                 
170                 fhPhiNeutral  = new TH2F
171                 ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T trigger}",
172                 nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
173                 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
174                 fhPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
175                 
176                 fhEtaNeutral  = new TH2F
177                 ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T trigger}",
178                  nptbins,ptmin,ptmax,netabins,etamin,etamax); 
179                 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
180                 fhEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
181                 
182                 fhDeltaPhiNeutral  = new TH2F
183                 ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
184                 nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
185                 fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
186                 fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
187                 
188                 fhDeltaPhiNeutralPt  = new TH2F
189                 ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
190                  nptbins,ptmin,ptmax,200,0,6.4); 
191                 fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
192                 fhDeltaPhiNeutralPt->SetXTitle("p_{T trigger} (GeV/c)");
193                 
194                 fhDeltaEtaNeutral  = new TH2F
195                 ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
196                  nptbins,ptmin,ptmax,200,-2,2); 
197                 fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
198                 fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
199                 
200                 fhPtImbalanceNeutral  = 
201                 new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
202                                  nptbins,ptmin,ptmax,1000,0.,1.2); 
203                 fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
204                 fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
205                 
206                 outputContainer->Add(fhPhiNeutral) ;  
207                 outputContainer->Add(fhEtaNeutral) ;   
208                 outputContainer->Add(fhDeltaPhiNeutral) ; 
209                 outputContainer->Add(fhDeltaEtaNeutral) ; 
210                 outputContainer->Add(fhPtImbalanceNeutral) ;
211                 
212                 
213                 //Keep neutral meson selection histograms if requiered
214                 //Setting done in AliNeutralMesonSelection
215                 if(GetNeutralMesonSelection()){
216                         TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
217                         if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
218                                 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
219                 }
220                 
221         }//Correlation with neutral hadrons
222         
223         return outputContainer;
224 }
225
226  //____________________________________________________________________________
227 void AliAnaParticleHadronCorrelation::InitParameters()
228 {
229  
230   //Initialize the parameters of the analysis.
231   SetInputAODName("photons");
232   SetPtCutRange(2,300);
233   fDeltaPhiMinCut = 1.5 ;
234   fDeltaPhiMaxCut = 4.5 ;
235   fSelectIsolated = kFALSE;
236 }
237
238 //__________________________________________________________________
239 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
240 {
241
242   //Print some relevant parameters set for the analysis
243   if(! opt)
244     return;
245   
246   Info("*** Print *** ", "%s %s", GetName(), GetTitle() ) ;
247   printf("pT Hadron       >    %2.2f\n", GetMinPt()) ; 
248   printf("pT Hadron       <    %2.2f\n", GetMaxPt()) ; 
249   printf("Phi trigger particle-Hadron      <     %3.2f\n", fDeltaPhiMaxCut) ; 
250   printf("Phi trigger particle-Hadron      >     %3.2f\n", fDeltaPhiMinCut) ;
251   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
252
253
254 //____________________________________________________________________________
255 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()  
256 {  
257         //Particle-Hadron Correlation Analysis, fill AODs
258         
259         if(!GetInputAODBranch())
260                 AliFatal(Form("ParticleHadronCorrelation::FillAOD: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data()));
261         
262         if(GetDebug() > 1){
263                 printf("Begin hadron correlation analysis, fill AODs \n");
264                 printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
265                 printf("In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
266                 printf("In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
267                 printf("In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
268         }
269         
270         //Loop on stored AOD particles, trigger
271         Int_t naod = GetInputAODBranch()->GetEntriesFast();
272         for(Int_t iaod = 0; iaod < naod ; iaod++){
273                 AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
274                 //Make correlation with charged hadrons
275                 if(GetReader()->IsCTSSwitchedOn() )
276                         MakeChargedCorrelation(particle, (TSeqCollection*)GetAODCTS(),kFALSE);
277                 
278                 //Make correlation with neutral pions
279                 //Trigger particle in PHOS, correlation with EMCAL
280                 if(particle->GetDetector()=="PHOS" && GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0)
281                         MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
282                 //Trigger particle in EMCAL, correlation with PHOS
283                 else if(particle->GetDetector()=="EMCAL" && GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0)
284                         MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
285                 //Trigger particle in CTS, correlation with PHOS, EMCAL and CTS
286                 else if(particle->GetDetector()=="CTS" ){
287                         if(GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0) 
288                                 MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
289                         if(GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0) 
290                                 MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
291                 }
292                 
293         }//Aod branch loop
294         
295         if(GetDebug() > 1) printf("End hadron correlation analysis, fill AODs \n");
296         
297 }
298
299 //____________________________________________________________________________
300 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()  
301 {  
302         //Particle-Hadron Correlation Analysis, fill histograms
303         
304         if(!GetInputAODBranch())
305                 AliFatal(Form("ParticleHadronCorrelation::FillHistos: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data()));
306         
307         if(GetDebug() > 1){
308                 printf("Begin hadron correlation analysis, fill histograms \n");
309                 printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
310         }
311         
312         //Loop on stored AOD particles
313         Int_t naod = GetInputAODBranch()->GetEntriesFast();
314         for(Int_t iaod = 0; iaod < naod ; iaod++){
315                 AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
316                 
317                 if(GetDebug() > 1){
318                         printf("Particle %d, In Track Refs  entries %d\n", iaod, (particle->GetRefTracks())->GetEntriesFast());
319                         printf("Particle %d, In Cluster Refs entries %d\n",iaod, (particle->GetRefClusters())->GetEntriesFast());
320                 }
321                 
322                 if(OnlyIsolated() && !particle->IsIsolated()) continue;
323                 
324                 //Make correlation with charged hadrons
325                 if((particle->GetRefTracks())->GetEntriesFast() > 0)
326                         MakeChargedCorrelation(particle, (TSeqCollection*) (particle->GetRefTracks()),kTRUE);
327                 
328                 //Make correlation with neutral pions
329                 if((particle->GetRefClusters())->GetEntriesFast() > 0)
330                         MakeNeutralCorrelation(particle,  (TSeqCollection*) (particle->GetRefClusters()), kTRUE);
331                 
332         }//Aod branch loop
333         
334         if(GetDebug() > 1) printf("End hadron correlation analysis, fill histograms \n");
335         
336 }
337
338 //____________________________________________________________________________
339 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TSeqCollection* pl, const Bool_t bFillHisto)
340 {  
341   // Charged Hadron Correlation Analysis
342   if(GetDebug() > 1)printf("Make trigger particle - charged hadron correlation \n");
343   
344   Double_t ptTrig  = aodParticle->Pt();
345   Double_t phiTrig = aodParticle->Phi();
346   Double_t pt   = -100.;
347   Double_t rat  = -100.; 
348   Double_t phi  = -100. ;
349   Double_t eta  = -100. ;
350   Double_t p[3];
351   
352   //Track loop, select tracks with good pt, phi and fill AODs or histograms
353   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
354     AliAODTrack * track = (AliAODTrack *) (pl->At(ipr)) ;
355     track->GetPxPyPz(p) ;
356     TLorentzVector mom(p[0],p[1],p[2],0);
357     pt    = mom.Pt();
358     eta  = mom.Eta();
359     phi  = mom.Phi() ;
360     if(phi<0) phi+=TMath::TwoPi();
361     rat   = pt/ptTrig ;
362
363     if(IsFidutialCutOn()){
364       Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
365       if(! in ) continue ;
366     }    
367
368     //Select only hadrons in pt range
369     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
370     
371     if(GetDebug() > 2)
372       printf("charged hadron: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max%2.2f, pT min %2.2f \n",
373              pt,phi,phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
374
375     if(bFillHisto){
376       // Fill Histograms
377       fhEtaCharged->Fill(ptTrig,eta);
378       fhPhiCharged->Fill(ptTrig,phi);
379       fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
380       fhDeltaPhiCharged->Fill(ptTrig,phiTrig-phi);
381       fhDeltaPhiChargedPt->Fill(pt,phiTrig-phi);
382       //Selection within angular range
383       if(((phiTrig-phi)> fDeltaPhiMinCut) && ((phiTrig-phi)<fDeltaPhiMaxCut) ){
384         if(GetDebug() > 2 ) printf("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f ",pt,phi,eta);
385         fhPtImbalanceCharged->Fill(ptTrig,rat);
386       } 
387     }
388     else{
389       //Fill AODs
390       aodParticle->AddTrack(track);
391       
392     }//aod particle loop
393   }// track loop
394
395 }  
396 //____________________________________________________________________________
397 void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle,TSeqCollection* pl, const Bool_t bFillHisto)  
398 {  
399   // Neutral Pion Correlation Analysis
400   if(GetDebug() > 1) printf("Make trigger particle - neutral hadron correlation \n");
401   
402   Double_t pt = -100.;
403   Double_t rat = -100.; 
404   Double_t phi = -100. ;
405   Double_t eta = -100. ;
406   Double_t ptTrig  = aodParticle->Pt();
407   Double_t phiTrig = aodParticle->Phi();
408   Double_t etaTrig = aodParticle->Eta();
409   
410   TLorentzVector gammai;
411   TLorentzVector gammaj;
412   
413   Double_t vertex[] = {0,0,0};
414
415   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
416   
417   //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
418   for(Int_t iclus = 0;iclus < pl->GetEntries() ; iclus ++ ){
419     AliAODCaloCluster * calo = (AliAODCaloCluster *) (pl->At(iclus)) ;
420     if(!bFillHisto){
421
422       //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
423       Int_t pdg=0;
424       if(!SelectCluster(calo, vertex, gammai, pdg)) continue ;
425       
426     if(GetDebug() > 2)
427       printf("neutral cluster: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max%2.2f, pT min %2.2f \n",
428              gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
429     
430       //2 gamma overlapped, found with PID
431       if(pdg == AliCaloPID::kPi0){ 
432
433         //Select only hadrons in pt range
434         if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) continue ;
435
436         aodParticle->AddCluster(calo);
437         if(GetDebug() > 2) printf("Correlated with selected pi0 (pid): pt %f, phi %f",gammai.Pt(),gammai.Phi());
438
439       }// pdg = 111
440
441       //Make invariant mass analysis
442       else if(pdg == AliCaloPID::kPhoton){      
443         //Search the photon companion in case it comes from  a Pi0 decay
444         //Apply several cuts to select the good pair;
445         for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
446           AliAODCaloCluster * calo2 = (AliAODCaloCluster *) (pl->At(jclus)) ;
447           
448           //Cluster selection, not charged with photon or pi0 id and in fidutial cut
449           Int_t pdgj=0;
450           if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
451
452           if(pdgj == AliCaloPID::kPhoton ){
453             
454             if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) continue ;
455             
456             //Select good pair (aperture and invariant mass)
457             if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
458
459               if(GetDebug() > 2 ) printf("Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f",
460                                          (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
461               Int_t labels[]={calo->GetLabel(0),calo2->GetLabel(0)};
462               Float_t pid[]={0,0,0,0,0,0,1,0,0,0,0,0,0};//Pi0 weight 1
463               Float_t pos[]={(gammai+gammaj).X(), (gammai+gammaj).Y(), (gammai+gammaj).Z()};
464
465               AliAODCaloCluster *caloCluster =  new AliAODCaloCluster(0,2,labels,(gammai+gammaj).E(), pos, pid,calo->GetType(),0);
466               aodParticle->AddCluster(caloCluster);
467             }//Pair selected
468           }//if pair of gammas
469         }//2nd loop
470       }// if pdg = 22
471     }// Fill AODs
472     else{ //Fill histograms
473       
474       calo->GetMomentum(gammai,vertex);//Assume that come from vertex in straight line
475       pt  = gammai.Pt();
476       
477       if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
478       
479       rat = pt/ptTrig ;
480       phi = gammai.Phi() ;
481       eta = gammai.Eta() ;
482
483       if(GetDebug() > 2 ) printf("Neutral Hadron Correlation: Histograms selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta);
484       
485       fhEtaNeutral->Fill(ptTrig,eta);
486       fhPhiNeutral->Fill(ptTrig,phi);
487       fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
488       fhDeltaPhiNeutral->Fill(ptTrig,phiTrig-phi);
489       fhDeltaPhiNeutralPt->Fill(pt,phiTrig-phi);
490       //Selection within angular range
491       if(((phiTrig-phi)> fDeltaPhiMinCut) && ((phiTrig-phi)<fDeltaPhiMaxCut) ){
492         if(GetDebug() > 2 ) printf("Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f ",pt,phi,eta);
493         fhPtImbalanceNeutral->Fill(ptTrig,rat);
494       }    
495     }//Fill histograms
496   }//1st loop
497   
498 }
499
500 //____________________________________________________________________________
501 Bool_t  AliAnaParticleHadronCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) const {
502    //Select cluster depending on its pid and acceptance selections
503    
504    //Skip matched clusters with tracks
505   if(calo->GetNTracksMatched() > 0) {
506   return kFALSE;} 
507    
508    //Check PID
509    calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
510    pdg = AliCaloPID::kPhoton;   
511    if(IsCaloPIDOn()){
512      //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
513      //or redo PID, recommended option for EMCal.
514      TString detector = "";
515      if(calo->IsPHOSCluster()) detector= "PHOS";
516      else if(calo->IsEMCALCluster()) detector= "EMCAL";
517
518      if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
519        pdg = GetCaloPID()->GetPdg(detector,calo->PID(),mom.E());//PID with weights
520      else
521        pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
522      
523      if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
524      
525      //If it does not pass pid, skip
526      if(pdg != AliCaloPID::kPhoton || pdg != AliCaloPID::kPi0) {
527        return kFALSE ;
528      }
529    }
530    
531    //Check acceptance selection
532    if(IsFidutialCutOn()){
533      Bool_t in = kFALSE;
534      if(calo->IsPHOSCluster())
535        in = GetFidutialCut()->IsInFidutialCut(mom,"PHOS") ;
536      else if(calo->IsEMCALCluster())
537        in = GetFidutialCut()->IsInFidutialCut(mom,"EMCAL") ;
538      if(! in ) { return kFALSE ;}
539    }
540    
541    if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
542    
543    return kTRUE;
544    
545  }