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()->GetEntriesFast());
 250     printf("In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
 251     printf("In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
 252     printf("In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
 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()->GetEntriesFast() > 0)
 266       MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
 267     //Trigger particle in EMCAL, correlation with PHOS
 268     else if(particle->GetDetector()=="EMCAL" && GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 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()->GetEntriesFast() > 0)
 273         MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
 274       if(GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 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()->GetEntriesFast());
 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())->GetEntriesFast());
 300       printf("Particle %d, In Cluster Refs entries %d\n",iaod, (particle->GetRefClusters())->GetEntriesFast());
 301     }
 302
 303     //Make correlation with charged hadrons
 304     if((particle->GetRefTracks())->GetEntriesFast() > 0)
 305       MakeChargedCorrelation(particle, (TSeqCollection*) (particle->GetRefTracks()),kTRUE);
 306
 307     //Make correlation with neutral pions
 308     if((particle->GetRefClusters())->GetEntriesFast() > 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  }