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  }