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