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  }