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