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