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