-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-/* $Id: $ */
-
-//_________________________________________________________________________
-// Class for the analysis of particle - hadron correlations
-// Particle (for example direct gamma) must be found in a previous analysis
-//-- Author: Gustavo Conesa (LNF-INFN)
-
-// Modified by Yaxian Mao:
-// 1. add the UE subtraction for corrlation study
-// 2. change the correlation variable
-// 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)
-// 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)
-// 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06)
-//////////////////////////////////////////////////////////////////////////////
-
-
-// --- ROOT system ---
-//#include "TClonesArray.h"
-#include "TClass.h"
-#include "TMath.h"
-#include "TH3D.h"
-
-//---- ANALYSIS system ----
-#include "AliNeutralMesonSelection.h"
-#include "AliAnaParticleHadronCorrelation.h"
-#include "AliCaloTrackReader.h"
-#include "AliCaloPID.h"
-#include "AliAODPWG4ParticleCorrelation.h"
-#include "AliFiducialCut.h"
-#include "AliVTrack.h"
-#include "AliVCluster.h"
-#include "AliMCAnalysisUtils.h"
-#include "TParticle.h"
-#include "AliStack.h"
-#include "AliAODMCParticle.h"
-#include "AliMixedEvent.h"
-
-ClassImp(AliAnaParticleHadronCorrelation)
-
-
-//____________________________________________________________________________
- AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
- AliAnaPartCorrBaseClass(),
- fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
- fMakeSeveralUE(0), fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.),
- fPi0AODBranchName(""),fNeutralCorr(0), fPi0Trigger(0),
- // fMultiBin(0),fNZvertBin(0),fNrpBin(0),fZvtxCut(0.),
- // fUseSelectEvent(kFALSE),
- // fhNclustersNtracks(0), //fhVertex(0),
- fhPtLeading(0),fhPhiLeading(0),fhEtaLeading(0),
- fhDeltaPhiDeltaEtaCharged(0),
- fhPhiCharged(0), fhEtaCharged(0),
- fhDeltaPhiCharged(0),
- fhDeltaEtaCharged(0),
- fhDeltaPhiChargedPt(0),
- fhDeltaPhiUeChargedPt(0),
- fhPtImbalanceCharged(0),
- fhPtImbalanceUeCharged(0),
- fhPtImbalancePosCharged(0),fhPtImbalanceNegCharged(0),
- fhPtHbpCharged(0), fhPtHbpUeCharged(0),
- fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),
- fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),
- fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0),
- fhPoutPtTrigPtAssoc(0), fhUePoutPtTrigPtAssoc(0),
- fhPtTrigCharged(0),
- fhTrigDeltaPhiDeltaEtaCharged(0x0), fhTrigCorr(0x0),fhTrigUeCorr(0x0),
- fhDeltaPhiDeltaEtaNeutral(0),
- fhPhiNeutral(0), fhEtaNeutral(0),
- fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),
- fhDeltaPhiNeutralPt(0),fhDeltaPhiUeNeutralPt(0),
- fhPtImbalanceNeutral(0),fhPtImbalanceUeNeutral(0),
- fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),
- fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),
- fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),
- fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0),
- fhPtPi0DecayRatio(0),
- fhDeltaPhiDecayCharged(0),
- fhPtImbalanceDecayCharged(0),
- fhDeltaPhiDecayNeutral(0),
- fhPtImbalanceDecayNeutral(0)
-{
- //Default Ctor
-
- //Initialize parameters
- InitParameters();
-}
-
-//________________________________________________________________________
-TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
-{
-
- // Create histograms to be saved in output file and
- // store them in fOutputContainer
- TList * outputContainer = new TList() ;
- outputContainer->SetName("CorrelationHistos") ;
-
- Int_t nptbins = GetHistoPtBins();
- Int_t nphibins = GetHistoPhiBins();
- Int_t netabins = GetHistoEtaBins();
- Float_t ptmax = GetHistoPtMax();
- Float_t phimax = GetHistoPhiMax();
- Float_t etamax = GetHistoEtaMax();
- Float_t ptmin = GetHistoPtMin();
- Float_t phimin = GetHistoPhiMin();
- Float_t etamin = GetHistoEtaMin();
-
-
-// fhNclustersNtracks = new TH2F ("Multiplicity","Neutral cluster and charged track multiplicity",1000, 0., 1000.,1000, 0., 1000.);
-// fhNclustersNtracks->SetYTitle("# of tracks");
-// fhNclustersNtracks->SetXTitle("# of clusters");
-
- fhPtLeading = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
- fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
-
- fhPhiLeading = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
- fhPhiLeading->SetYTitle("#phi (rad)");
-
- fhEtaLeading = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
- fhEtaLeading->SetYTitle("#eta ");
- // outputContainer->Add(fhNclustersNtracks);
- outputContainer->Add(fhPtLeading);
- outputContainer->Add(fhPhiLeading);
- outputContainer->Add(fhEtaLeading);
-
- //Correlation with charged hadrons
- if(GetReader()->IsCTSSwitchedOn()) {
- fhDeltaPhiDeltaEtaCharged = new TH2F
- ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
- 140,-2.,5.,200,-2,2);
- fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
- fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
-
- fhPhiCharged = new TH2F
- ("PhiCharged","#phi_{h^{#pm}} vs p_{T #pm}",
- nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
- fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
-
- fhEtaCharged = new TH2F
- ("EtaCharged","#eta_{h^{#pm}} vs p_{T #pm}",
- nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
- fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
-
- fhDeltaPhiCharged = new TH2F
- ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiCharged->SetYTitle("#Delta #phi");
- fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
-
- fhDeltaPhiChargedPt = new TH2F
- ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
- fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
-
- fhDeltaPhiUeChargedPt = new TH2F
- ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
- fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
-
- fhDeltaEtaCharged = new TH2F
- ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
- nptbins,ptmin,ptmax,200,-2,2);
- fhDeltaEtaCharged->SetYTitle("#Delta #eta");
- fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
-
- fhPtImbalanceCharged =
- new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
- fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
-
- fhPtImbalanceUeCharged =
- new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
- fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
-
- fhPtImbalancePosCharged =
- new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");
- fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");
-
- fhPtImbalanceNegCharged =
- new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");
- fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");
-
- fhPtHbpCharged =
- new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
- fhPtHbpCharged->SetXTitle("p_{T trigger}");
-
- fhPtHbpUeCharged =
- new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
- fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
-
- fhPoutPtTrigPtAssoc =
- new TH3D("PoutPtTrigPtAssoc","Pout with charged hadrons",
- nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
- fhPoutPtTrigPtAssoc->SetZTitle("p_{out} (GeV/c)");
- fhPoutPtTrigPtAssoc->SetYTitle("p_{T associated} (GeV/c)");
- fhPoutPtTrigPtAssoc->SetXTitle("p_{T trigger} (GeV/c)");
-
- fhUePoutPtTrigPtAssoc =
- new TH3D("UePoutPtTrigPtAssoc"," UE Pout with charged hadrons",
- nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
- fhUePoutPtTrigPtAssoc->SetZTitle("p_{out} (GeV/c)");
- fhUePoutPtTrigPtAssoc->SetYTitle("p_{T associated} (GeV/c)");
- fhUePoutPtTrigPtAssoc->SetXTitle("p_{T trigger} (GeV/c)");
-
- fhPtTrigCharged =
- new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",
- nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
- fhPtTrigCharged->SetYTitle("p_{T asso} (GeV/c)");
- fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
-
- outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
- outputContainer->Add(fhPhiCharged) ;
- outputContainer->Add(fhEtaCharged) ;
- outputContainer->Add(fhDeltaPhiCharged) ;
- outputContainer->Add(fhDeltaEtaCharged) ;
- outputContainer->Add(fhDeltaPhiChargedPt) ;
- outputContainer->Add(fhDeltaPhiUeChargedPt) ;
- outputContainer->Add(fhPtImbalanceCharged) ;
- outputContainer->Add(fhPtImbalancePosCharged) ;
- outputContainer->Add(fhPtImbalanceNegCharged) ;
- outputContainer->Add(fhPtImbalanceUeCharged) ;
- outputContainer->Add(fhPtHbpCharged) ;
- outputContainer->Add(fhPtHbpUeCharged) ;
- outputContainer->Add(fhPoutPtTrigPtAssoc) ;
- outputContainer->Add(fhUePoutPtTrigPtAssoc) ;
- outputContainer->Add(fhPtTrigCharged) ;
-
- if(DoEventSelect()){
- Int_t nMultiBins = GetMultiBin();
- fhTrigDeltaPhiDeltaEtaCharged = new TH3D*[nMultiBins] ;
- fhTrigCorr = new TH2F*[nMultiBins];
- fhTrigUeCorr = new TH2F*[nMultiBins];
- for(Int_t im=0; im<nMultiBins; im++){
- fhTrigDeltaPhiDeltaEtaCharged[im] = new TH3D
- (Form("fhTrigDeltaPhiDeltaEtaCharged_%d",im),Form("fhTrigDeltaPhiDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5., 200,-2,2);
- fhTrigDeltaPhiDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
- fhTrigDeltaPhiDeltaEtaCharged[im]->SetYTitle("#Delta #phi");
- fhTrigDeltaPhiDeltaEtaCharged[im]->SetZTitle("#Delta #eta");
- fhTrigCorr[im] = new TH2F
- (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
- fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
- fhTrigCorr[im]->SetXTitle("p_{T trigger}");
- fhTrigUeCorr[im] = new TH2F
- (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
- fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
- fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");
-
- outputContainer->Add(fhTrigDeltaPhiDeltaEtaCharged[im]) ;
- outputContainer->Add(fhTrigCorr[im]);
- outputContainer->Add(fhTrigUeCorr[im]);
-
- }
- }
-
- if(fPi0Trigger){
- fhPtPi0DecayRatio = new TH2F
- ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
- nptbins,ptmin,ptmax, 100,0.,2.);
- fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
- fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
-
- fhDeltaPhiDecayCharged = new TH2F
- ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
- fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
-
- fhPtImbalanceDecayCharged =
- new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");
- fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");
-
- outputContainer->Add(fhPtPi0DecayRatio) ;
- outputContainer->Add(fhDeltaPhiDecayCharged) ;
- outputContainer->Add(fhPtImbalanceDecayCharged) ;
- }
-
-
- if(fMakeSeveralUE){
- fhDeltaPhiUeLeftCharged = new TH2F
- ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
- fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
- outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
-
- fhDeltaPhiUeRightCharged = new TH2F
- ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
- fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
- outputContainer->Add(fhDeltaPhiUeRightCharged) ;
-
- fhPtImbalanceUeLeftCharged =
- new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
- fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
-
- fhPtImbalanceUeRightCharged =
- new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
- fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtImbalanceUeRightCharged) ;
-
- fhPtHbpUeLeftCharged =
- new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
- fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtHbpUeLeftCharged) ;
-
- fhPtHbpUeRightCharged =
- new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
- fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtHbpUeRightCharged) ;
-
- }
- } //Correlation with charged hadrons
-
- //Correlation with neutral hadrons
- if(fNeutralCorr){
-
- fhDeltaPhiDeltaEtaNeutral = new TH2F
- ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
- 140,-2.,5.,200,-2,2);
- fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
- fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
-
- fhPhiNeutral = new TH2F
- ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #pi^{0}}",
- nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
- fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
-
- fhEtaNeutral = new TH2F
- ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #pi^{0}}",
- nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
- fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
-
- fhDeltaPhiNeutral = new TH2F
- ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
- nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
- fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
-
- fhDeltaPhiNeutralPt = new TH2F
- ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
- fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
-
- fhDeltaPhiUeNeutralPt = new TH2F
- ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
- fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
-
- fhDeltaEtaNeutral = new TH2F
- ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
- nptbins,ptmin,ptmax,200,-2,2);
- fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
- fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
-
- fhPtImbalanceNeutral =
- new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
- fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
-
- fhPtImbalanceUeNeutral =
- new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
- fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
-
- fhPtHbpNeutral =
- new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
- fhPtHbpNeutral->SetXTitle("p_{T trigger}");
-
- fhPtHbpUeNeutral =
- new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
- fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
-
-
- outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
- outputContainer->Add(fhPhiNeutral) ;
- outputContainer->Add(fhEtaNeutral) ;
- outputContainer->Add(fhDeltaPhiNeutral) ;
- outputContainer->Add(fhDeltaPhiNeutralPt) ;
- outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
- outputContainer->Add(fhDeltaEtaNeutral) ;
- outputContainer->Add(fhPtImbalanceNeutral) ;
- outputContainer->Add(fhPtImbalanceUeNeutral) ;
- outputContainer->Add(fhPtHbpNeutral) ;
- outputContainer->Add(fhPtHbpUeNeutral) ;
-
- if(fPi0Trigger){
- fhDeltaPhiDecayNeutral = new TH2F
- ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
- fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
-
- fhPtImbalanceDecayNeutral =
- new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");
- fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");
-
- outputContainer->Add(fhDeltaPhiDecayNeutral) ;
- outputContainer->Add(fhPtImbalanceDecayNeutral) ;
- }
-
-
- if(fMakeSeveralUE){
- fhDeltaPhiUeLeftNeutral = new TH2F
- ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
- fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
- outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
-
- fhDeltaPhiUeRightNeutral = new TH2F
- ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
- nptbins,ptmin,ptmax,140,-2.,5.);
- fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
- fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
- outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
-
- fhPtImbalanceUeLeftNeutral =
- new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
- nptbins,ptmin,ptmax,140,0.,2.);
- fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
- fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
-
- fhPtImbalanceUeRightNeutral =
- new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
- nptbins,ptmin,ptmax,200,0.,2.);
- fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
- fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
-
- fhPtHbpUeLeftNeutral =
- new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
- fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtHbpUeLeftNeutral) ;
-
- fhPtHbpUeRightNeutral =
- new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
- nptbins,ptmin,ptmax,200,0.,10.);
- fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
- fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
- outputContainer->Add(fhPtHbpUeRightNeutral) ;
-
- }
-
- //Keep neutral meson selection histograms if requiered
- //Setting done in AliNeutralMesonSelection
- if(GetNeutralMesonSelection()){
- TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
- if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
- for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
- delete nmsHistos;
- }
-
- }//Correlation with neutral hadrons
-
- return outputContainer;
-
-}
-
-//____________________________________________________________________________
-void AliAnaParticleHadronCorrelation::InitParameters()
-{
-
- //Initialize the parameters of the analysis.
- SetInputAODName("PWG4Particle");
- SetAODObjArrayName("Hadrons");
- AddToHistogramsName("AnaHadronCorr_");
-
- SetPtCutRange(0.,300);
- fDeltaPhiMinCut = 1.5 ;
- fDeltaPhiMaxCut = 4.5 ;
- fSelectIsolated = kFALSE;
- fMakeSeveralUE = kFALSE;
- fUeDeltaPhiMinCut = 1. ;
- fUeDeltaPhiMaxCut = 1.5 ;
- fNeutralCorr = kFALSE ;
- fPi0Trigger = kFALSE ;
-
-}
-
-//__________________________________________________________________
-void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
-{
-
- //Print some relevant parameters set for the analysis
- if(! opt)
- return;
-
- printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
- AliAnaPartCorrBaseClass::Print(" ");
-
- printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ;
- printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;
- printf("Isolated Trigger? %d\n", fSelectIsolated) ;
- printf("Phi trigger particle-UeHadron < %3.2f\n", fUeDeltaPhiMaxCut) ;
- printf("Phi trigger particle-UeHadron > %3.2f\n", fUeDeltaPhiMinCut) ;
- printf("Several UE? %d\n", fMakeSeveralUE) ;
- printf("Name of AOD Pi0 Branch %s \n",fPi0AODBranchName.Data());
- printf("Do Decay-hadron correlation ? %d\n", fPi0Trigger) ;
-
-
-}
-
-//__________________________________________________________________
-TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
-{
- //Save parameters used for analysis
- TString parList ; //this will be list of parameters used for this analysis.
- const Int_t buffersize = 255;
- char onePar[buffersize] ;
-
- snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Phi trigger particle-Hadron < %3.2f ", fDeltaPhiMaxCut) ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Phi trigger particle-Hadron > %3.2f ", fDeltaPhiMinCut) ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Isolated Trigger? %d\n", fSelectIsolated) ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Phi trigger particle-UeHadron < %3.2f ", fUeDeltaPhiMaxCut) ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Phi trigger particle-UeHadron > %3.2f ", fUeDeltaPhiMinCut) ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Several UE? %d\n", fMakeSeveralUE) ;
- parList+=onePar ;
- snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ",fPi0AODBranchName.Data());
- parList+=onePar ;
- snprintf(onePar,buffersize,"Do Decay-hadron correlation ? %d", fPi0Trigger) ;
- parList+=onePar ;
-
- //Get parameters set in base class.
- parList += GetBaseParametersList() ;
-
- //Get parameters set in PID class.
- //parList += GetCaloPID()->GetPIDParametersList() ;
-
- //Get parameters set in FiducialCut class (not available yet)
- //parlist += GetFidCut()->GetFidCutParametersList()
-
- return new TObjString(parList) ;
-
-}
-
-
-//____________________________________________________________________________
-void AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
-{
- //Particle-Hadron Correlation Analysis, fill AODs
-
- if(!GetInputAODBranch()){
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
- abort();
- }
-
- if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
- abort();
- }
-
- if(GetDebug() > 1){
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetPHOSClusters()->GetEntriesFast());
- }
-
- //Loop on stored AOD particles, trigger
- Double_t ptTrig = 0.;
- Int_t trigIndex = -1;
- Int_t naod = GetInputAODBranch()->GetEntriesFast();
- //fhNclustersNtracks->Fill(naod, GetCTSTracks()->GetEntriesFast());
- for(Int_t iaod = 0; iaod < naod ; iaod++){
- AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
- //find the leading particles with highest momentum
- if (particle->Pt()>ptTrig) {
- ptTrig = particle->Pt() ;
- trigIndex = iaod ;
- }
- }//Aod branch loop
-
- //Do correlation with leading particle
- if(trigIndex!=-1){
-
- AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
- //Make correlation with charged hadrons
- if(GetReader()->IsCTSSwitchedOn() )
- MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);
-
- TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
- if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
- MakeNeutralCorrelation(particle, pi0list,kFALSE);
-
- }//Correlate leading
-
- if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
-
-}
-
-//____________________________________________________________________________
-void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
-{
- //Particle-Hadron Correlation Analysis, fill histograms
-
- if(!GetInputAODBranch()){
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
- abort();
- }
-
- if(GetDebug() > 1){
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
- printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
- }
-
- //Get the vertex and check it is not too large in z
- Double_t v[3] = {0,0,0}; //vertex ;
- GetReader()->GetVertex(v);
- if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
-
- //Loop on stored AOD particles, find leading
- Int_t naod = GetInputAODBranch()->GetEntriesFast();
- // if(naod!=0)fhVertex->Fill(v[0],v[1],v[2]);
- Double_t ptTrig = 0.;
- Int_t trigIndex = -1 ;
- for(Int_t iaod = 0; iaod < naod ; iaod++){ //loop on input trigger AOD file
- AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
- //vertex cut in case of mixing
- if (GetMixedEvent()) {
- Int_t evt=-1;
- Int_t id =-1;
- if (particle->GetCaloLabel(0) >= 0 ){
- id=particle->GetCaloLabel(0);
- if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
- }
- else if(particle->GetTrackLabel(0) >= 0 ){
- id=particle->GetTrackLabel(0);
- if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
- }
- else continue;
-
- if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
- return ;
- }
-
- //check if the particle is isolated or if we want to take the isolation into account
- if(OnlyIsolated() && !particle->IsIsolated()) continue;
- //check if inside the vertex cut
- //find the leading particles with highest momentum
- if (particle->Pt()>ptTrig) {
- ptTrig = particle->Pt() ;
- trigIndex = iaod ;
- }
- }//finish searching for leading trigger particle
- if(trigIndex!=-1){ //using trigger partilce to do correlations
- AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
-
- if (GetMixedEvent()) {
- Int_t evt=-1;
- Int_t id = 0;
- if (particle->GetCaloLabel(0) >= 0 ){
- id=particle->GetCaloLabel(0);
- if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
- }
- else if(particle->GetTrackLabel(0) >= 0 ){
- id=particle->GetTrackLabel(0);
- if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
- }
- else return;
-
- if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
- return ;
- }
-
- //Fill leading particle histogram
- fhPtLeading->Fill(particle->Pt());
- Float_t phi = particle->Phi();
- if(phi<0)phi+=TMath::TwoPi();
- fhPhiLeading->Fill(particle->Pt(), phi);
- fhEtaLeading->Fill(particle->Pt(), particle->Eta());
-
- //Make correlation with charged hadrons
- if(GetReader()->IsCTSSwitchedOn() )
- MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);
-
- TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
- if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
- MakeNeutralCorrelation(particle, pi0list,kTRUE);
-
- }//Aod branch loop
-
- if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
-
-}
-
-//____________________________________________________________________________
-void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)
-{
- // Charged Hadron Correlation Analysis
- if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
-
- Int_t evtIndex11 = -1 ; //cluster trigger or pi0 trigger
- Int_t evtIndex12 = -1 ; // pi0 trigger
- Int_t evtIndex13 = -1 ; // charged trigger
- Int_t indexPhoton1 = -1 ;
- Int_t indexPhoton2 = -1 ;
- Double_t v[3] = {0,0,0}; //vertex ;
- GetReader()->GetVertex(v);
- if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
-
- Int_t nTracks = GetCTSTracks()->GetEntriesFast() ;
-
- if (GetMixedEvent()) {
- evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
- evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
- evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
- }
-
- Double_t ptTrig = aodParticle->Pt();
- Double_t pxTrig = aodParticle->Px();
- Double_t pyTrig = aodParticle->Py();
-
- Double_t phiTrig = aodParticle->Phi();
- Double_t etaTrig = aodParticle->Eta();
-
- Double_t pt = -100.;
- Double_t px = -100.;
- Double_t py = -100.;
- Double_t rat = -100.;
- Double_t xE = -100.;
- Double_t cosi = -100.;
- Double_t phi = -100. ;
- Double_t eta = -100. ;
- TVector3 p3;
-
- TObjArray * reftracks = 0x0;
- Int_t nrefs = 0;
-
- Double_t ptDecay1 = 0. ;
- Double_t pxDecay1 = 0. ;
- Double_t pyDecay1 = 0. ;
- Double_t phiDecay1 = 0. ;
- Double_t ptDecay2 = 0. ;
- Double_t pxDecay2 = 0. ;
- Double_t pyDecay2 = 0. ;
- Double_t phiDecay2 = 0. ;
-
- Double_t ratDecay1 = -100.;
- Double_t ratDecay2 = -100.;
- Float_t deltaphi = -100. ;
- Float_t deltaphiDecay1 = -100. ;
- Float_t deltaphiDecay2 = -100. ;
- TObjArray * clusters = 0x0 ;
- TLorentzVector photonMom ;
-
- if(fPi0Trigger){
- indexPhoton1 = aodParticle->GetCaloLabel (0);
- indexPhoton2 = aodParticle->GetCaloLabel (1);
- if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
-
- if(indexPhoton1!=-1 && indexPhoton2!=-1){
- if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
- else clusters = GetPHOSClusters() ;
- for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
- AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
- photon->GetMomentum(photonMom,GetVertex(0)) ;
- if(photon->GetID()==indexPhoton1) {
- ptDecay1 = photonMom.Pt();
- pxDecay1 = photonMom.Px();
- pyDecay1 = photonMom.Py();
- phiDecay1 = photonMom.Phi();
- if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);
- }
- if(photon->GetID()==indexPhoton2) {
- ptDecay2 = photonMom.Pt();
- pxDecay2 = photonMom.Px();
- pyDecay2 = photonMom.Py();
- phiDecay2 = photonMom.Phi();
- if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);
- }
- if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
- } //cluster loop
- } //index of decay photons found
-
- } //make decay-hadron correlation
-
- //Track loop, select tracks with good pt, phi and fill AODs or histograms
- //Int_t currentIndex = -1 ;
- for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
- AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
-
- //check if inside the vertex cut
- //printf("charge = %d\n", track->Charge());
- Int_t evtIndex2 = 0 ;
- if (GetMixedEvent()) {
- evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
- if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
- continue ;
- //vertex cut
- if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
- return ;
- // if(currentIndex == evtIndex2) // tracks from different event
- // continue ;
- // currentIndex = evtIndex2 ;
- }
- Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
- p3.SetXYZ(mom[0],mom[1],mom[2]);
- pt = p3.Pt();
- px = p3.Px();
- py = p3.Py();
- eta = p3.Eta();
- phi = p3.Phi() ;
- if(phi < 0) phi+=TMath::TwoPi();
-
- //Select only hadrons in pt range
- if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
- //remove trigger itself for correlation when use charged triggers
- if(track->GetID()==aodParticle->GetTrackLabel(0) && pt==ptTrig && phi==phiTrig && eta==etaTrig)
- continue ;
- if(IsFiducialCutOn()){
- Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
- if(! in ) continue ;
- }
- //jumped out this event if near side associated partile pt larger than trigger
- if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;
- rat = pt/ptTrig ;
- xE = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
- if(xE <0.)xE =-xE;
- cosi = TMath::Log(1/xE);
- // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
- // printf("phi = %f \n", phi);
-
- if(fPi0Trigger){
- if(indexPhoton1!=-1 && indexPhoton2!=-1){
- if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
- if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
- deltaphiDecay1 = phiDecay1-phi;
- deltaphiDecay2 = phiDecay2-phi;
- if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();
- if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();
- if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();
- if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();
- }
- } //do decay-hadron correlation
-
- //Selection within angular range
- deltaphi = phiTrig-phi;
- if(deltaphi< -TMath::PiOver2()) deltaphi+=TMath::TwoPi();
- if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
-
- Double_t pout = pt*TMath::Sin(deltaphi) ;
-
- if(GetDebug() > 2)
- printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
- pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());
-
- if(bFillHisto){
- // Fill Histograms
- fhEtaCharged->Fill(pt,eta);
- fhPhiCharged->Fill(pt,phi);
- fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
- fhDeltaPhiCharged->Fill(ptTrig, deltaphi);
- fhDeltaPhiDeltaEtaCharged->Fill(deltaphi,aodParticle->Eta()-eta);
-
- if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
- //fill different multiplicity histogram
- if(DoEventSelect()){
- for(Int_t im=0; im<GetMultiBin(); im++){
- if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))fhTrigDeltaPhiDeltaEtaCharged[im]->Fill(ptTrig,deltaphi,aodParticle->Eta()-eta);
- }
- }
- //delta phi cut for correlation
- if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
- fhDeltaPhiChargedPt->Fill(pt,deltaphi);
- fhPtImbalanceCharged->Fill(ptTrig,xE);
- fhPtHbpCharged->Fill(ptTrig,cosi);
- fhPoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
- fhPtTrigCharged->Fill(ptTrig, pt) ;
- if(track->Charge()>0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;
- else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;
- //fill different multiplicity histogram
- if(DoEventSelect()){
- for(Int_t im=0; im<GetMultiBin(); im++){
- if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
- fhTrigCorr[im]->Fill(ptTrig,xE);
- }
- } //multiplicity events selection
- } //delta phi cut for correlation
- else { //UE study
- fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);
- fhUePoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
- fhPtImbalanceUeCharged->Fill(ptTrig,xE);
- fhPtHbpUeCharged->Fill(ptTrig,cosi);
- if(DoEventSelect()){
- for(Int_t im=0; im<GetMultiBin(); im++){
- if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
- fhTrigUeCorr[im]->Fill(ptTrig,xE);
- }
- } //multiplicity events selection
-
- } //UE study
-
- if(fPi0Trigger){
- if(indexPhoton1!=-1 && indexPhoton2!=-1){
- fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaphiDecay1);
- fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaphiDecay2);
- if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
- if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
- fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1);
- if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
- fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);
- if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
- } //index of decay photons found
- } //make decay-hadron correlation
-
- //several UE calculation
- if(fMakeSeveralUE){
- if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
- fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);
- fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
- fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
- }
- if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
- fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);
- fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
- fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
-
- }
- } //several UE calculation
-
- } //Fill histogram
- else{
- nrefs++;
- if(nrefs==1){
- reftracks = new TObjArray(0);
- reftracks->SetName(GetAODObjArrayName()+"Tracks");
- reftracks->SetOwner(kFALSE);
- }
- reftracks->Add(track);
- }//aod particle loop
- }// track loop
-
- //Fill AOD with reference tracks, if not filling histograms
- if(!bFillHisto && reftracks) {
- aodParticle->AddObjArray(reftracks);
- }
-
- //delete reftracks;
-
-}
-
-
-
-//____________________________________________________________________________
-//void AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector)
-//{
-// // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed
-// if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");
-//
-// if(!NewOutputAOD()){
-// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");
-// abort();
-// }
-//
-// Double_t phiTrig = aodParticle->Phi();
-// Int_t tag = 0;
-// TLorentzVector gammai;
-// TLorentzVector gammaj;
-//
-// //Get vertex for photon momentum calculation
-//
-// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
-// {
-// for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {
-// if (!GetMixedEvent())
-// GetReader()->GetVertex(GetVertex(iev));
-// else
-// GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev));
-// }
-// }
-// Double_t vertex2[] = {0.0,0.0,0.0} ; //vertex of second input aod
-// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
-// {
-// if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
-// }
-//
-// //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
-// //Int_t iEvent= GetReader()->GetEventNumber() ;
-// Int_t nclus = pl->GetEntriesFast();
-// for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){
-// AliVCluster * calo = (AliVCluster *) (pl->At(iclus)) ;
-//
-// Int_t evtIndex1 = 0 ;
-// if (GetMixedEvent()) {
-// evtIndex1=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
-// }
-//
-//
-// //Input from second AOD?
-// Int_t inputi = 0;
-// if (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= iclus)
-// inputi = 1 ;
-// else if(aodParticle->GetDetector() == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= iclus)
-// inputi = 1;
-//
-// //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
-// //FIXME
-// Int_t pdg=0;
-// //if (inputi == 0 && !SelectCluster(calo, GetVertex(evtIndex1), gammai, pdg))
-// continue ;
-// //MEFIX
-// else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg))
-// continue ;
-//
-// if(GetDebug() > 2)
-// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral cluster in %s: pt %f, phi %f, phi trigger %f. Cuts: delta phi min %2.2f, max %2.2f, pT min %2.2f \n",
-// detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
-//
-// //2 gamma overlapped, found with PID
-// if(pdg == AliCaloPID::kPi0){
-//
-// //Select only hadrons in pt range
-// if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt())
-// continue ;
-//
-// //Selection within angular range
-// Float_t phi = gammai.Phi();
-// if(phi < 0) phi+=TMath::TwoPi();
-// //Float_t deltaphi = TMath::Abs(phiTrig-phi);
-// //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
-//
-// AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);
-// //pi0.SetLabel(calo->GetLabel());
-// pi0.SetPdg(AliCaloPID::kPi0);
-// pi0.SetDetector(detector);
-//
-// if(IsDataMC()){
-// pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(),inputi));
-// if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());
-// }//Work with stack also
-// //Set the indeces of the original caloclusters
-// pi0.SetCaloLabel(calo->GetID(),-1);
-// AddAODParticle(pi0);
-//
-// if(GetDebug() > 2)
-// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());
-//
-// }// pdg = 111
-//
-// //Make invariant mass analysis
-// else if(pdg == AliCaloPID::kPhoton){
-// //Search the photon companion in case it comes from a Pi0 decay
-// //Apply several cuts to select the good pair;
-// for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
-// AliVCluster * calo2 = (AliVCluster *) (pl->At(jclus)) ;
-// Int_t evtIndex2 = 0 ;
-// if (GetMixedEvent()) {
-// evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
-// }
-// if (GetMixedEvent() && (evtIndex1 == evtIndex2))
-// continue ;
-//
-// //Input from second AOD?
-// Int_t inputj = 0;
-// if (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= jclus)
-// inputj = 1;
-// else if(aodParticle->GetDetector() == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= jclus)
-// inputj = 1;
-//
-// //Cluster selection, not charged with photon or pi0 id and in fiducial cut
-// Int_t pdgj=0;
-// //FIXME
-// //if (inputj == 0 && !SelectCluster(calo2, GetVertex(evtIndex2), gammaj, pdgj))
-// continue ;
-// //MEFIX
-//
-// else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))
-// continue ;
-// //FIXME
-// //if(!SelectCluster(calo2,GetVertex(evtIndex2), gammaj, pdgj))
-// //MEFIX
-// continue ;
-//
-// if(pdgj == AliCaloPID::kPhoton ){
-//
-// if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt())
-// continue ;
-//
-// //Selection within angular range
-// Float_t phi = (gammai+gammaj).Phi();
-// if(phi < 0) phi+=TMath::TwoPi();
-// //Float_t deltaphi = TMath::Abs(phiTrig-phi);
-// //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
-//
-// //Select good pair (aperture and invariant mass)
-// if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
-//
-// if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
-// (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
-//
-// TLorentzVector pi0mom = gammai+gammaj;
-// AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);
-// //pi0.SetLabel(calo->GetLabel());
-// pi0.SetPdg(AliCaloPID::kPi0);
-// pi0.SetDetector(detector);
-// if(IsDataMC()){
-// //Check origin of the candidates
-//
-// Int_t label1 = calo->GetLabel();
-// Int_t label2 = calo2->GetLabel();
-// Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);
-// Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);
-//
-// if(GetDebug() > 0)
-// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
-// if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
-//
-// //Check if pi0 mother is the same
-// if(GetReader()->ReadStack()){
-// TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
-// label1 = mother1->GetFirstMother();
-// //mother1 = GetMCStack()->Particle(label1);//pi0
-//
-// TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
-// label2 = mother2->GetFirstMother();
-// //mother2 = GetMCStack()->Particle(label2);//pi0
-// }
-// else if(GetReader()->ReadAODMCParticles()){
-// AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree
-// label1 = mother1->GetMother();
-// //mother1 = GetMCStack()->Particle(label1);//pi0
-// AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree
-// label2 = mother2->GetMother();
-// //mother2 = GetMCStack()->Particle(label2);//pi0
-// }
-//
-// //printf("mother1 %d, mother2 %d\n",label1,label2);
-// if(label1 == label2)
-// GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
-// }
-// }//Work with mc information also
-// pi0.SetTag(tag);
-// //Set the indeces of the original caloclusters
-// pi0.SetCaloLabel(calo->GetID(), calo2->GetID());
-// AddAODParticle(pi0);
-//
-//
-// }//Pair selected
-// }//if pair of gammas
-// }//2nd loop
-// }// if pdg = 22
-// }//1st loop
-//
-// if(GetDebug() > 1)
-// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());
-//}
-
-//____________________________________________________________________________
-void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, TObjArray* pi0list, const Bool_t bFillHisto)
-{
- // Neutral Pion Correlation Analysis
- if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());
-
- Int_t evtIndex11 = 0 ;
- Int_t evtIndex12 = 0 ;
- if (GetMixedEvent()) {
- evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
- evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
- }
-
- Double_t pt = -100.;
- Double_t px = -100.;
- Double_t py = -100.;
- Double_t rat = -100.;
- Double_t phi = -100.;
- Double_t eta = -100.;
- Double_t xE = -100.;
- Double_t cosi = -100.;
-
- Double_t ptTrig = aodParticle->Pt();
- Double_t phiTrig = aodParticle->Phi();
- Double_t etaTrig = aodParticle->Eta();
- Double_t pxTrig = aodParticle->Px();
- Double_t pyTrig = aodParticle->Py();
-
-
- Int_t indexPhoton1 = -1 ;
- Int_t indexPhoton2 = -1 ;
- Double_t ptDecay1 = 0. ;
- Double_t pxDecay1 = 0. ;
- Double_t pyDecay1 = 0. ;
- Double_t phiDecay1 = 0. ;
- Double_t ptDecay2 = 0. ;
- Double_t pxDecay2 = 0. ;
- Double_t pyDecay2 = 0. ;
- Double_t phiDecay2 = 0. ;
-
- Double_t ratDecay1 = -100.;
- Double_t ratDecay2 = -100.;
- Float_t deltaphi = -100. ;
- Float_t deltaphiDecay1 = -100. ;
- Float_t deltaphiDecay2 = -100. ;
- TObjArray * clusters = 0x0 ;
- TLorentzVector photonMom ;
- if(fPi0Trigger){
- indexPhoton1 = aodParticle->GetCaloLabel (0);
- indexPhoton2 = aodParticle->GetCaloLabel (1);
- if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
-
- if(indexPhoton1!=-1 && indexPhoton2!=-1){
- if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
- else clusters = GetPHOSClusters() ;
- for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
- AliVCluster * photon = (AliVCluster*) (clusters->At(iclus));
- photon->GetMomentum(photonMom,GetVertex(0)) ;
- if(photon->GetID()==indexPhoton1) {
- ptDecay1 = photonMom.Pt();
- pxDecay1 = photonMom.Px();
- pyDecay1 = photonMom.Py();
- phiDecay1 = photonMom.Phi();
- }
- if(photon->GetID()==indexPhoton2) {
- ptDecay2 = photonMom.Pt();
- pxDecay2 = photonMom.Px();
- pyDecay2 = photonMom.Py();
- phiDecay2 = photonMom.Phi();
- }
- if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
- } //photonAOD loop
- } //index of decay photons found
- if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);
- } //make decay-hadron correlation
-
- TObjArray * refpi0 =0x0;
- Int_t nrefs = 0;
-
- //Loop on stored AOD pi0
- Int_t naod = pi0list->GetEntriesFast();
- if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - aod branch entries %d\n", naod);
- for(Int_t iaod = 0; iaod < naod ; iaod++){
- AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (pi0list->At(iaod));
-
- Int_t evtIndex2 = 0 ;
- Int_t evtIndex3 = 0 ;
- if (GetMixedEvent()) {
- evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
- evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
-
- if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
- continue ;
- }
-
- //Int_t pdg = pi0->GetPdg();
- //if(pdg != AliCaloPID::kPi0) continue;
-
- pt = pi0->Pt();
- px = pi0->Px();
- py = pi0->Py();
- if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
- //jumped out this event if near side associated partile pt larger than trigger
- if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;
-
- //Selection within angular range
- phi = pi0->Phi();
- //Float_t deltaphi = TMath::Abs(phiTrig-phi);
- //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
-
- if(bFillHisto){
-
- deltaphi = phiTrig-phi;
- if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();
- if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
-
- rat = pt/ptTrig ;
- phi = pi0->Phi() ;
- eta = pi0->Eta() ;
- xE = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
- if(xE <0.)xE =-xE;
- cosi = TMath::Log(1/xE);
-
- if(fPi0Trigger){
- if(indexPhoton1!=-1 && indexPhoton2!=-1){
- if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
- if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
- deltaphiDecay1 = phiDecay1-phi;
- deltaphiDecay2 = phiDecay2-phi;
- if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();
- if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();
- if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();
- if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();
- fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaphiDecay1);
- fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaphiDecay2);
- if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
- if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
- fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1);
- if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
- fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);
- if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
- }
- } //do decay-hadron correlation
-
- fhEtaNeutral->Fill(pt,eta);
- fhPhiNeutral->Fill(pt,phi);
- fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
- fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);
- fhDeltaPhiDeltaEtaNeutral->Fill(deltaphi,etaTrig-eta);
-
- //delta phi cut for correlation
- if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
- fhDeltaPhiNeutralPt->Fill(pt,deltaphi);
- fhPtImbalanceNeutral->Fill(ptTrig,rat);
- fhPtHbpNeutral->Fill(ptTrig,cosi);
- }
- else {
- fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);
- fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
- fhPtHbpUeNeutral->Fill(ptTrig,cosi);
- }
- //several UE calculation
- if(fMakeSeveralUE){
- if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
- fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);
- fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
- fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
- }
- if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
- fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);
- fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
- fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
- }
- } //several UE calculation
- }
- else{
- nrefs++;
- if(nrefs==1){
- refpi0 = new TObjArray(0);
- refpi0->SetName(GetAODObjArrayName()+"Pi0s");
- refpi0->SetOwner(kFALSE);
- }
- refpi0->Add(pi0);
- }//put references in trigger AOD
-
- //if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
-
- }//loop
-}
-
-
-//____________________________________________________________________________
-//Bool_t AliAnaParticleHadronCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {
-// //Select cluster depending on its pid and acceptance selections
-//
-// //Skip matched clusters with tracks
-// if(IsTrackMatched(calo)) return kFALSE;
-//
-// TString detector = "";
-// if (calo->IsPHOS()) detector= "PHOS";
-// else if(calo->IsEMCAL()) detector= "EMCAL";
-//
-// //Check PID
-// calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
-// pdg = AliCaloPID::kPhoton;
-// if(IsCaloPIDOn()){
-// //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
-// //or redo PID, recommended option for EMCal.
-//
-// if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
-// pdg = GetCaloPID()->GetPdg(detector,calo->GetPID(),mom.E());//PID with weights
-// else
-// pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
-//
-// if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
-//
-// //If it does not pass pid, skip
-// if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {
-// return kFALSE ;
-// }
-// }//PID on
-//
-// //Check acceptance selection
-// if(IsFiducialCutOn()){
-// Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,detector) ;
-// if(! in ) return kFALSE ;
-// }
-//
-// if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
-//
-// return kTRUE;
-//
-//}
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes is hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+/* $Id: $ */\r
+\r
+//_________________________________________________________________________\r
+// Class for the analysis of particle - hadron correlations\r
+// Particle (for example direct gamma) must be found in a previous analysis \r
+//-- Author: Gustavo Conesa (LNF-INFN) \r
+\r
+// Modified by Yaxian Mao:\r
+// 1. add the UE subtraction for corrlation study\r
+// 2. change the correlation variable\r
+// 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)\r
+// 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)\r
+// 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06) \r
+//////////////////////////////////////////////////////////////////////////////\r
+\r
+\r
+// --- ROOT system ---\r
+//#include "TClonesArray.h"\r
+#include "TClass.h"\r
+#include "TMath.h"\r
+#include "TH3D.h"\r
+\r
+//---- ANALYSIS system ----\r
+#include "AliNeutralMesonSelection.h" \r
+#include "AliAnaParticleHadronCorrelation.h" \r
+#include "AliCaloTrackReader.h"\r
+#include "AliCaloPID.h"\r
+#include "AliAODPWG4ParticleCorrelation.h"\r
+#include "AliFiducialCut.h"\r
+#include "AliVTrack.h"\r
+#include "AliVCluster.h"\r
+#include "AliMCAnalysisUtils.h"\r
+#include "TParticle.h"\r
+#include "AliStack.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliMixedEvent.h"\r
+\r
+ClassImp(AliAnaParticleHadronCorrelation)\r
+\r
+\r
+//____________________________________________________________________________\r
+ AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(): \r
+ AliAnaPartCorrBaseClass(),\r
+ fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),\r
+ fMakeSeveralUE(0), \r
+ fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.), \r
+ fPi0AODBranchName(""),fNeutralCorr(0), fPi0Trigger(0),\r
+ // fMultiBin(0),fNZvertBin(0),fNrpBin(0),fZvtxCut(0.),\r
+ // fUseSelectEvent(kFALSE), \r
+ // fhNclustersNtracks(0), //fhVertex(0),\r
+ fhPtLeading(0),fhPhiLeading(0),fhEtaLeading(0), \r
+ fhDeltaPhiDeltaEtaCharged(0),\r
+ fhPhiCharged(0), fhEtaCharged(0), \r
+ fhDeltaPhiCharged(0), \r
+ fhDeltaEtaCharged(0), \r
+ fhDeltaPhiChargedPt(0), \r
+ fhDeltaPhiUeChargedPt(0), \r
+ fhPtImbalanceCharged(0), \r
+ fhPtImbalanceUeCharged(0),\r
+ fhPtImbalancePosCharged(0),fhPtImbalanceNegCharged(0),\r
+ fhPtHbpCharged(0), fhPtHbpUeCharged(0),\r
+ fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),\r
+ fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),\r
+ fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0), \r
+ fhPtTrigPout(0), fhPtAssocDeltaPhi(0),\r
+ fhPtTrigCharged(0),\r
+ fhTrigDeltaPhiCharged(0x0), fhTrigDeltaEtaCharged(0x0),fhTrigCorr(0x0),fhTrigUeCorr(0x0),\r
+ fhDeltaPhiDeltaEtaNeutral(0), \r
+ fhPhiNeutral(0), fhEtaNeutral(0), \r
+ fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),\r
+ fhDeltaPhiNeutralPt(0),fhDeltaPhiUeNeutralPt(0), \r
+ fhPtImbalanceNeutral(0),fhPtImbalanceUeNeutral(0),\r
+ fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),\r
+ fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),\r
+ fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),\r
+ fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0),\r
+ fhPtPi0DecayRatio(0),\r
+ fhDeltaPhiDecayCharged(0),\r
+ fhPtImbalanceDecayCharged(0), \r
+ fhDeltaPhiDecayNeutral(0),\r
+ fhPtImbalanceDecayNeutral(0)\r
+{\r
+ //Default Ctor\r
+ \r
+ //Initialize parameters\r
+ InitParameters();\r
+}\r
+\r
+//________________________________________________________________________\r
+TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()\r
+{ \r
+ \r
+ // Create histograms to be saved in output file and \r
+ // store them in fOutputContainer\r
+ TList * outputContainer = new TList() ; \r
+ outputContainer->SetName("CorrelationHistos") ; \r
+ \r
+ Int_t nptbins = GetHistoPtBins();\r
+ Int_t nphibins = GetHistoPhiBins();\r
+ Int_t netabins = GetHistoEtaBins();\r
+ Float_t ptmax = GetHistoPtMax();\r
+ Float_t phimax = GetHistoPhiMax();\r
+ Float_t etamax = GetHistoEtaMax();\r
+ Float_t ptmin = GetHistoPtMin();\r
+ Float_t phimin = GetHistoPhiMin();\r
+ Float_t etamin = GetHistoEtaMin(); \r
+\r
+ \r
+// fhNclustersNtracks = new TH2F ("Multiplicity","Neutral cluster and charged track multiplicity",1000, 0., 1000.,1000, 0., 1000.); \r
+// fhNclustersNtracks->SetYTitle("# of tracks");\r
+// fhNclustersNtracks->SetXTitle("# of clusters");\r
+\r
+ fhPtLeading = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax); \r
+ fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");\r
+ \r
+ fhPhiLeading = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax); \r
+ fhPhiLeading->SetYTitle("#phi (rad)");\r
+ \r
+ fhEtaLeading = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax); \r
+ fhEtaLeading->SetYTitle("#eta "); \r
+ // outputContainer->Add(fhNclustersNtracks);\r
+ outputContainer->Add(fhPtLeading);\r
+ outputContainer->Add(fhPhiLeading);\r
+ outputContainer->Add(fhEtaLeading);\r
+ \r
+ //Correlation with charged hadrons\r
+ if(GetReader()->IsCTSSwitchedOn()) {\r
+ fhDeltaPhiDeltaEtaCharged = new TH2F\r
+ ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",\r
+ 140,-2.,5.,200,-2,2); \r
+ fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");\r
+ fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta"); \r
+ \r
+ fhPhiCharged = new TH2F\r
+ ("PhiCharged","#phi_{h^{#pm}} vs p_{T #pm}",\r
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
+ fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");\r
+ fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");\r
+ \r
+ fhEtaCharged = new TH2F\r
+ ("EtaCharged","#eta_{h^{#pm}} vs p_{T #pm}",\r
+ nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
+ fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");\r
+ fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");\r
+ \r
+ fhDeltaPhiCharged = new TH2F\r
+ ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.); \r
+ fhDeltaPhiCharged->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");\r
+ \r
+ fhDeltaPhiChargedPt = new TH2F\r
+ ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.);\r
+ fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
+ \r
+ fhDeltaPhiUeChargedPt = new TH2F\r
+ ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.);\r
+ fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
+ \r
+ fhDeltaEtaCharged = new TH2F\r
+ ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,-2,2); \r
+ fhDeltaEtaCharged->SetYTitle("#Delta #eta");\r
+ fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");\r
+ \r
+ fhPtImbalanceCharged = \r
+ new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");\r
+ fhPtImbalanceCharged->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtImbalanceUeCharged = \r
+ new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");\r
+ fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtImbalancePosCharged = \r
+ new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");\r
+ fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtImbalanceNegCharged = \r
+ new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");\r
+ fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtHbpCharged = \r
+ new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpCharged->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpCharged->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtHbpUeCharged = \r
+ new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpUeCharged->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtTrigPout = \r
+ new TH2F("PtTrigPout","Pout with triggers",\r
+ nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax); \r
+ fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");\r
+ fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)"); \r
+ \r
+ fhPtAssocDeltaPhi = \r
+ new TH2F("fhPtAssocDeltaPhi"," charged hadrons vs. delta phi",\r
+ nptbins,ptmin,ptmax,140,-2.,5.); \r
+ fhPtAssocDeltaPhi->SetXTitle("p_{T h^{#pm}} (GeV/c)"); \r
+ fhPtAssocDeltaPhi->SetYTitle("#Delta #phi (GeV/c)"); \r
+ \r
+ fhPtTrigCharged = \r
+ new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",\r
+ nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); \r
+ fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");\r
+ fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)"); \r
+ \r
+ outputContainer->Add(fhDeltaPhiDeltaEtaCharged);\r
+ outputContainer->Add(fhPhiCharged) ;\r
+ outputContainer->Add(fhEtaCharged) ;\r
+ outputContainer->Add(fhDeltaPhiCharged) ; \r
+ outputContainer->Add(fhDeltaEtaCharged) ;\r
+ outputContainer->Add(fhDeltaPhiChargedPt) ;\r
+ outputContainer->Add(fhDeltaPhiUeChargedPt) ;\r
+ outputContainer->Add(fhPtImbalanceCharged) ;\r
+ outputContainer->Add(fhPtImbalancePosCharged) ;\r
+ outputContainer->Add(fhPtImbalanceNegCharged) ;\r
+ outputContainer->Add(fhPtImbalanceUeCharged) ;\r
+ outputContainer->Add(fhPtHbpCharged) ;\r
+ outputContainer->Add(fhPtHbpUeCharged) ;\r
+ outputContainer->Add(fhPtTrigPout) ;\r
+ outputContainer->Add(fhPtAssocDeltaPhi) ;\r
+ outputContainer->Add(fhPtTrigCharged) ;\r
+ \r
+ if(DoEventSelect()){ \r
+ Int_t nMultiBins = GetMultiBin();\r
+ fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;\r
+ fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;\r
+ fhTrigCorr = new TH2F*[nMultiBins];\r
+ fhTrigUeCorr = new TH2F*[nMultiBins];\r
+ for(Int_t im=0; im<nMultiBins; im++){\r
+ fhTrigDeltaPhiCharged[im] = new TH2F \r
+ (Form("fhTrigDeltaPhiCharged_%d",im),Form("fhTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.); \r
+ fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");\r
+ fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");\r
+ fhTrigDeltaEtaCharged[im] = new TH2F \r
+ (Form("fhTrigDeltaEtaCharged_%d",im),Form("fhTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2); \r
+ fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");\r
+ fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");\r
+ fhTrigCorr[im] = new TH2F\r
+ (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");\r
+ fhTrigCorr[im]->SetXTitle("p_{T trigger}");\r
+ fhTrigUeCorr[im] = new TH2F\r
+ (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");\r
+ fhTrigUeCorr[im]->SetXTitle("p_{T trigger}"); \r
+\r
+ outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;\r
+ outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;\r
+ outputContainer->Add(fhTrigCorr[im]);\r
+ outputContainer->Add(fhTrigUeCorr[im]);\r
+ \r
+ }\r
+ }\r
+ \r
+ if(fPi0Trigger){\r
+ fhPtPi0DecayRatio = new TH2F\r
+ ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay", \r
+ nptbins,ptmin,ptmax, 100,0.,2.); \r
+ fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");\r
+ fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");\r
+ \r
+ fhDeltaPhiDecayCharged = new TH2F\r
+ ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.); \r
+ fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");\r
+ \r
+ fhPtImbalanceDecayCharged = \r
+ new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");\r
+ fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");\r
+ \r
+ outputContainer->Add(fhPtPi0DecayRatio) ; \r
+ outputContainer->Add(fhDeltaPhiDecayCharged) ; \r
+ outputContainer->Add(fhPtImbalanceDecayCharged) ;\r
+ } \r
+ \r
+ \r
+ if(fMakeSeveralUE){ \r
+ fhDeltaPhiUeLeftCharged = new TH2F\r
+ ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",\r
+ nptbins,ptmin,ptmax,140,-2.,5.);\r
+ fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
+ outputContainer->Add(fhDeltaPhiUeLeftCharged) ;\r
+ \r
+ fhDeltaPhiUeRightCharged = new TH2F\r
+ ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",\r
+ nptbins,ptmin,ptmax,140,-2.,5.);\r
+ fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
+ outputContainer->Add(fhDeltaPhiUeRightCharged) ;\r
+ \r
+ fhPtImbalanceUeLeftCharged = \r
+ new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");\r
+ fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtImbalanceUeLeftCharged) ;\r
+ \r
+ fhPtImbalanceUeRightCharged = \r
+ new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");\r
+ fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtImbalanceUeRightCharged) ;\r
+ \r
+ fhPtHbpUeLeftCharged = \r
+ new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtHbpUeLeftCharged) ;\r
+ \r
+ fhPtHbpUeRightCharged = \r
+ new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtHbpUeRightCharged) ;\r
+ \r
+ } \r
+ } //Correlation with charged hadrons\r
+ \r
+ //Correlation with neutral hadrons\r
+ if(fNeutralCorr){\r
+ \r
+ fhDeltaPhiDeltaEtaNeutral = new TH2F\r
+ ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",\r
+ 140,-2.,5.,200,-2,2); \r
+ fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");\r
+ fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta"); \r
+ \r
+ fhPhiNeutral = new TH2F\r
+ ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #pi^{0}}",\r
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
+ fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");\r
+ fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
+ \r
+ fhEtaNeutral = new TH2F\r
+ ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #pi^{0}}",\r
+ nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
+ fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");\r
+ fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
+ \r
+ fhDeltaPhiNeutral = new TH2F\r
+ ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",\r
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
+ fhDeltaPhiNeutral->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");\r
+ \r
+ fhDeltaPhiNeutralPt = new TH2F\r
+ ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.); \r
+ fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");\r
+ \r
+ fhDeltaPhiUeNeutralPt = new TH2F\r
+ ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.); \r
+ fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");\r
+ \r
+ fhDeltaEtaNeutral = new TH2F\r
+ ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,-2,2); \r
+ fhDeltaEtaNeutral->SetYTitle("#Delta #eta");\r
+ fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");\r
+ \r
+ fhPtImbalanceNeutral = \r
+ new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");\r
+ fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtImbalanceUeNeutral = \r
+ new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");\r
+ fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtHbpNeutral = \r
+ new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpNeutral->SetXTitle("p_{T trigger}");\r
+ \r
+ fhPtHbpUeNeutral = \r
+ new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");\r
+ \r
+ \r
+ outputContainer->Add(fhDeltaPhiDeltaEtaNeutral); \r
+ outputContainer->Add(fhPhiNeutral) ; \r
+ outputContainer->Add(fhEtaNeutral) ; \r
+ outputContainer->Add(fhDeltaPhiNeutral) ; \r
+ outputContainer->Add(fhDeltaPhiNeutralPt) ; \r
+ outputContainer->Add(fhDeltaPhiUeNeutralPt) ; \r
+ outputContainer->Add(fhDeltaEtaNeutral) ; \r
+ outputContainer->Add(fhPtImbalanceNeutral) ;\r
+ outputContainer->Add(fhPtImbalanceUeNeutral) ; \r
+ outputContainer->Add(fhPtHbpNeutral) ;\r
+ outputContainer->Add(fhPtHbpUeNeutral) ; \r
+ \r
+ if(fPi0Trigger){\r
+ fhDeltaPhiDecayNeutral = new TH2F\r
+ ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",\r
+ nptbins,ptmin,ptmax,140,-2.,5.); \r
+ fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");\r
+ \r
+ fhPtImbalanceDecayNeutral = \r
+ new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");\r
+ fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");\r
+ \r
+ outputContainer->Add(fhDeltaPhiDecayNeutral) ; \r
+ outputContainer->Add(fhPtImbalanceDecayNeutral) ;\r
+ }\r
+ \r
+ \r
+ if(fMakeSeveralUE){ \r
+ fhDeltaPhiUeLeftNeutral = new TH2F\r
+ ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",\r
+ nptbins,ptmin,ptmax,140,-2.,5.);\r
+ fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");\r
+ outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;\r
+ \r
+ fhDeltaPhiUeRightNeutral = new TH2F\r
+ ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",\r
+ nptbins,ptmin,ptmax,140,-2.,5.);\r
+ fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");\r
+ fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");\r
+ outputContainer->Add(fhDeltaPhiUeRightNeutral) ;\r
+ \r
+ fhPtImbalanceUeLeftNeutral = \r
+ new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",\r
+ nptbins,ptmin,ptmax,140,0.,2.); \r
+ fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");\r
+ fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;\r
+ \r
+ fhPtImbalanceUeRightNeutral = \r
+ new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,2.); \r
+ fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");\r
+ fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtImbalanceUeRightNeutral) ;\r
+ \r
+ fhPtHbpUeLeftNeutral = \r
+ new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtHbpUeLeftNeutral) ;\r
+ \r
+ fhPtHbpUeRightNeutral = \r
+ new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",\r
+ nptbins,ptmin,ptmax,200,0.,10.); \r
+ fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");\r
+ fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");\r
+ outputContainer->Add(fhPtHbpUeRightNeutral) ;\r
+ \r
+ } \r
+ \r
+ //Keep neutral meson selection histograms if requiered\r
+ //Setting done in AliNeutralMesonSelection\r
+ if(GetNeutralMesonSelection()){\r
+ TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;\r
+ if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())\r
+ for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;\r
+ delete nmsHistos;\r
+ }\r
+ \r
+ }//Correlation with neutral hadrons\r
+ \r
+ return outputContainer;\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaParticleHadronCorrelation::InitParameters()\r
+{\r
+ \r
+ //Initialize the parameters of the analysis.\r
+ SetInputAODName("PWG4Particle");\r
+ SetAODObjArrayName("Hadrons"); \r
+ AddToHistogramsName("AnaHadronCorr_");\r
+\r
+ SetPtCutRange(0.,300);\r
+ fDeltaPhiMinCut = 1.5 ;\r
+ fDeltaPhiMaxCut = 4.5 ;\r
+ fSelectIsolated = kFALSE;\r
+ fMakeSeveralUE = kFALSE;\r
+ fUeDeltaPhiMinCut = 1. ;\r
+ fUeDeltaPhiMaxCut = 1.5 ;\r
+ fNeutralCorr = kFALSE ;\r
+ fPi0Trigger = kFALSE ;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const\r
+{\r
+\r
+ //Print some relevant parameters set for the analysis\r
+ if(! opt)\r
+ return;\r
+ \r
+ printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
+ AliAnaPartCorrBaseClass::Print(" ");\r
+\r
+ printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ; \r
+ printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;\r
+ printf("Isolated Trigger? %d\n", fSelectIsolated) ;\r
+ printf("Phi trigger particle-UeHadron < %3.2f\n", fUeDeltaPhiMaxCut) ; \r
+ printf("Phi trigger particle-UeHadron > %3.2f\n", fUeDeltaPhiMinCut) ;\r
+ printf("Several UE? %d\n", fMakeSeveralUE) ;\r
+ printf("Name of AOD Pi0 Branch %s \n",fPi0AODBranchName.Data());\r
+ printf("Do Decay-hadron correlation ? %d\n", fPi0Trigger) ;\r
+\r
+ \r
+} \r
+\r
+//__________________________________________________________________\r
+TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()\r
+{\r
+ //Save parameters used for analysis\r
+ TString parList ; //this will be list of parameters used for this analysis.\r
+ const Int_t buffersize = 255;\r
+ char onePar[buffersize] ;\r
+ \r
+ snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"Phi trigger particle-Hadron < %3.2f ", fDeltaPhiMaxCut) ; \r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Phi trigger particle-Hadron > %3.2f ", fDeltaPhiMinCut) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Isolated Trigger? %d\n", fSelectIsolated) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Phi trigger particle-UeHadron < %3.2f ", fUeDeltaPhiMaxCut) ; \r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Phi trigger particle-UeHadron > %3.2f ", fUeDeltaPhiMinCut) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Several UE? %d\n", fMakeSeveralUE) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ",fPi0AODBranchName.Data());\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"Do Decay-hadron correlation ? %d", fPi0Trigger) ;\r
+ parList+=onePar ;\r
+\r
+ //Get parameters set in base class.\r
+ parList += GetBaseParametersList() ;\r
+ \r
+ //Get parameters set in PID class.\r
+ //parList += GetCaloPID()->GetPIDParametersList() ;\r
+ \r
+ //Get parameters set in FiducialCut class (not available yet)\r
+ //parlist += GetFidCut()->GetFidCutParametersList() \r
+ \r
+ return new TObjString(parList) ; \r
+\r
+} \r
+\r
+\r
+//____________________________________________________________________________\r
+void AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() \r
+{ \r
+ //Particle-Hadron Correlation Analysis, fill AODs\r
+ \r
+ if(!GetInputAODBranch()){\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());\r
+ abort();\r
+ }\r
+ \r
+ if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());\r
+ abort();\r
+ }\r
+ \r
+ if(GetDebug() > 1){\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetPHOSClusters()->GetEntriesFast());\r
+ }\r
+ \r
+ //Loop on stored AOD particles, trigger\r
+ Double_t ptTrig = 0.;\r
+ Int_t trigIndex = -1;\r
+ Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
+ //fhNclustersNtracks->Fill(naod, GetCTSTracks()->GetEntriesFast());\r
+ for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+ AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
+ //find the leading particles with highest momentum\r
+ if (particle->Pt()>ptTrig) {\r
+ ptTrig = particle->Pt() ;\r
+ trigIndex = iaod ;\r
+ }\r
+ }//Aod branch loop\r
+ \r
+ //Do correlation with leading particle\r
+ if(trigIndex!=-1){\r
+ \r
+ AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));\r
+ //Make correlation with charged hadrons\r
+ if(GetReader()->IsCTSSwitchedOn() )\r
+ MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);\r
+ \r
+ TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists\r
+ if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)\r
+ MakeNeutralCorrelation(particle, pi0list,kFALSE);\r
+ \r
+ }//Correlate leading\r
+ \r
+ if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() \r
+{ \r
+ //Particle-Hadron Correlation Analysis, fill histograms\r
+ \r
+ if(!GetInputAODBranch()){\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());\r
+ abort();\r
+ }\r
+ \r
+ if(GetDebug() > 1){\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");\r
+ printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());\r
+ }\r
+ \r
+ //Get the vertex and check it is not too large in z\r
+ Double_t v[3] = {0,0,0}; //vertex ;\r
+ GetReader()->GetVertex(v);\r
+ if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ; \r
+ \r
+ //Loop on stored AOD particles, find leading\r
+ Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
+ // if(naod!=0)fhVertex->Fill(v[0],v[1],v[2]);\r
+ Double_t ptTrig = 0.;\r
+ Int_t trigIndex = -1 ;\r
+ for(Int_t iaod = 0; iaod < naod ; iaod++){ //loop on input trigger AOD file \r
+ AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
+ //vertex cut in case of mixing\r
+ if (GetMixedEvent()) {\r
+ Int_t evt=-1;\r
+ Int_t id =-1;\r
+ if (particle->GetCaloLabel(0) >= 0 ){\r
+ id=particle->GetCaloLabel(0); \r
+ if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;\r
+ }\r
+ else if(particle->GetTrackLabel(0) >= 0 ){\r
+ id=particle->GetTrackLabel(0);\r
+ if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;\r
+ }\r
+ else continue;\r
+ \r
+ if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut()) \r
+ return ;\r
+ }\r
+\r
+ //check if the particle is isolated or if we want to take the isolation into account\r
+ if(OnlyIsolated() && !particle->IsIsolated()) continue;\r
+ //check if inside the vertex cut\r
+ //find the leading particles with highest momentum\r
+ if (particle->Pt()>ptTrig) {\r
+ ptTrig = particle->Pt() ;\r
+ trigIndex = iaod ;\r
+ }\r
+ }//finish searching for leading trigger particle\r
+ if(trigIndex!=-1){ //using trigger partilce to do correlations\r
+ AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));\r
+ \r
+ if (GetMixedEvent()) {\r
+ Int_t evt=-1;\r
+ Int_t id = 0;\r
+ if (particle->GetCaloLabel(0) >= 0 ){\r
+ id=particle->GetCaloLabel(0); \r
+ if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;\r
+ }\r
+ else if(particle->GetTrackLabel(0) >= 0 ){\r
+ id=particle->GetTrackLabel(0);\r
+ if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;\r
+ }\r
+ else return;\r
+ \r
+ if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut()) \r
+ return ;\r
+ }\r
+ \r
+ //Fill leading particle histogram \r
+ fhPtLeading->Fill(particle->Pt());\r
+ Float_t phi = particle->Phi();\r
+ if(phi<0)phi+=TMath::TwoPi();\r
+ fhPhiLeading->Fill(particle->Pt(), phi);\r
+ fhEtaLeading->Fill(particle->Pt(), particle->Eta());\r
+ \r
+ //Make correlation with charged hadrons\r
+ if(GetReader()->IsCTSSwitchedOn() )\r
+ MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);\r
+ \r
+ TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists\r
+ if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)\r
+ MakeNeutralCorrelation(particle, pi0list,kTRUE);\r
+ \r
+ }//Aod branch loop\r
+ \r
+ if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)\r
+{ \r
+ // Charged Hadron Correlation Analysis\r
+ if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");\r
+ \r
+ Int_t evtIndex11 = -1 ; //cluster trigger or pi0 trigger \r
+ Int_t evtIndex12 = -1 ; // pi0 trigger\r
+ Int_t evtIndex13 = -1 ; // charged trigger\r
+ Int_t indexPhoton1 = -1 ;\r
+ Int_t indexPhoton2 = -1 ; \r
+ Double_t v[3] = {0,0,0}; //vertex ;\r
+ GetReader()->GetVertex(v);\r
+ if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ; \r
+\r
+ Int_t nTracks = GetCTSTracks()->GetEntriesFast() ;\r
+ \r
+ if (GetMixedEvent()) {\r
+ evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;\r
+ evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ; \r
+ evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;\r
+ }\r
+ \r
+ Double_t ptTrig = aodParticle->Pt();\r
+ Double_t pxTrig = aodParticle->Px();\r
+ Double_t pyTrig = aodParticle->Py();\r
+ \r
+ Double_t phiTrig = aodParticle->Phi();\r
+ Double_t etaTrig = aodParticle->Eta();\r
+ \r
+ Double_t pt = -100.;\r
+ Double_t px = -100.;\r
+ Double_t py = -100.;\r
+ Double_t rat = -100.; \r
+ Double_t xE = -100.; \r
+ Double_t cosi = -100.; \r
+ Double_t phi = -100. ;\r
+ Double_t eta = -100. ;\r
+ TVector3 p3; \r
+ \r
+ TObjArray * reftracks = 0x0;\r
+ Int_t nrefs = 0;\r
+ \r
+ Double_t ptDecay1 = 0. ;\r
+ Double_t pxDecay1 = 0. ;\r
+ Double_t pyDecay1 = 0. ;\r
+ Double_t phiDecay1 = 0. ;\r
+ Double_t ptDecay2 = 0. ;\r
+ Double_t pxDecay2 = 0. ;\r
+ Double_t pyDecay2 = 0. ;\r
+ Double_t phiDecay2 = 0. ;\r
+ \r
+ Double_t ratDecay1 = -100.; \r
+ Double_t ratDecay2 = -100.; \r
+ Float_t deltaphi = -100. ;\r
+ Float_t deltaphiDecay1 = -100. ;\r
+ Float_t deltaphiDecay2 = -100. ;\r
+ TObjArray * clusters = 0x0 ; \r
+ TLorentzVector photonMom ; \r
+ \r
+ if(fPi0Trigger){\r
+ indexPhoton1 = aodParticle->GetCaloLabel (0);\r
+ indexPhoton2 = aodParticle->GetCaloLabel (1);\r
+ if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);\r
+ \r
+ if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
+ if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;\r
+ else clusters = GetPHOSClusters() ;\r
+ for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){\r
+ AliVCluster * photon = (AliVCluster*) (clusters->At(iclus)); \r
+ photon->GetMomentum(photonMom,GetVertex(0)) ;\r
+ if(photon->GetID()==indexPhoton1) {\r
+ ptDecay1 = photonMom.Pt();\r
+ pxDecay1 = photonMom.Px();\r
+ pyDecay1 = photonMom.Py();\r
+ phiDecay1 = photonMom.Phi();\r
+ if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);\r
+ }\r
+ if(photon->GetID()==indexPhoton2) {\r
+ ptDecay2 = photonMom.Pt();\r
+ pxDecay2 = photonMom.Px();\r
+ pyDecay2 = photonMom.Py();\r
+ phiDecay2 = photonMom.Phi();\r
+ if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);\r
+ } \r
+ if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);\r
+ } //cluster loop \r
+ } //index of decay photons found\r
+ \r
+ } //make decay-hadron correlation\r
+ \r
+ //Track loop, select tracks with good pt, phi and fill AODs or histograms\r
+ //Int_t currentIndex = -1 ; \r
+ for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){\r
+ AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;\r
+\r
+ //check if inside the vertex cut\r
+ //printf("charge = %d\n", track->Charge());\r
+ Int_t evtIndex2 = 0 ; \r
+ if (GetMixedEvent()) {\r
+ evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;\r
+ if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events\r
+ continue ; \r
+ //vertex cut\r
+ if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) \r
+ return ;\r
+ // if(currentIndex == evtIndex2) // tracks from different event \r
+ // continue ;\r
+ // currentIndex = evtIndex2 ;\r
+ }\r
+ Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};\r
+ p3.SetXYZ(mom[0],mom[1],mom[2]);\r
+ pt = p3.Pt();\r
+ px = p3.Px();\r
+ py = p3.Py();\r
+ eta = p3.Eta();\r
+ phi = p3.Phi() ;\r
+ if(phi < 0) phi+=TMath::TwoPi();\r
+\r
+ //Select only hadrons in pt range\r
+ if(pt < GetMinPt() || pt > GetMaxPt()) continue ;\r
+ //remove trigger itself for correlation when use charged triggers \r
+ if(track->GetID()==aodParticle->GetTrackLabel(0) && pt==ptTrig && phi==phiTrig && eta==etaTrig) \r
+ continue ; \r
+ if(IsFiducialCutOn()){\r
+ Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;\r
+ if(! in ) continue ;\r
+ } \r
+ //jumped out this event if near side associated partile pt larger than trigger\r
+ if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;\r
+ rat = pt/ptTrig ;\r
+ xE = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);\r
+ if(xE <0.)xE =-xE;\r
+ cosi = TMath::Log(1/xE);\r
+ // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);\r
+ // printf("phi = %f \n", phi);\r
+ \r
+ if(fPi0Trigger){\r
+ if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
+ if(ptDecay1) ratDecay1 = pt/ptDecay1 ;\r
+ if(ptDecay2) ratDecay2 = pt/ptDecay2 ; \r
+ deltaphiDecay1 = phiDecay1-phi;\r
+ deltaphiDecay2 = phiDecay2-phi;\r
+ if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();\r
+ if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();\r
+ if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();\r
+ if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi(); \r
+ }\r
+ } //do decay-hadron correlation \r
+ \r
+ //Selection within angular range\r
+ deltaphi = phiTrig-phi;\r
+ if(deltaphi< -TMath::PiOver2()) deltaphi+=TMath::TwoPi();\r
+ if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();\r
+\r
+ Double_t pout = pt*TMath::Sin(deltaphi) ;\r
+ \r
+ if(GetDebug() > 2)\r
+ printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts: delta phi %2.2f < %2.2f < %2.2f, pT min %2.2f \n",\r
+ pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());\r
+ \r
+ if(bFillHisto){\r
+ // Fill Histograms\r
+ fhEtaCharged->Fill(pt,eta);\r
+ fhPhiCharged->Fill(pt,phi);\r
+ fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);\r
+ fhDeltaPhiCharged->Fill(ptTrig, deltaphi);\r
+ fhDeltaPhiDeltaEtaCharged->Fill(deltaphi,aodParticle->Eta()-eta);\r
+ fhPtAssocDeltaPhi->Fill(pt, deltaphi);\r
+ \r
+ if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);\r
+ //fill different multiplicity histogram\r
+ if(DoEventSelect()){\r
+ for(Int_t im=0; im<GetMultiBin(); im++){\r
+ if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1)){\r
+ fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaphi);\r
+ fhTrigDeltaEtaCharged[im]->Fill(ptTrig,aodParticle->Eta()-eta);\r
+ }\r
+ }\r
+ }\r
+ //delta phi cut for correlation\r
+ if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {\r
+ fhDeltaPhiChargedPt->Fill(pt,deltaphi);\r
+ fhPtImbalanceCharged->Fill(ptTrig,xE); \r
+ fhPtHbpCharged->Fill(ptTrig,cosi);\r
+ fhPtTrigPout->Fill(ptTrig, pout) ;\r
+ fhPtTrigCharged->Fill(ptTrig, pt) ;\r
+ if(track->Charge()>0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;\r
+ else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;\r
+ //fill different multiplicity histogram\r
+ if(DoEventSelect()){\r
+ for(Int_t im=0; im<GetMultiBin(); im++){\r
+ if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))\r
+ fhTrigCorr[im]->Fill(ptTrig,xE);\r
+ }\r
+ } //multiplicity events selection\r
+ } //delta phi cut for correlation\r
+ else if ((deltaphi > fUeDeltaPhiMinCut) && ( deltaphi < fUeDeltaPhiMaxCut)) { //UE study\r
+ fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);\r
+ //fhUePoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;\r
+ fhPtImbalanceUeCharged->Fill(ptTrig,xE);\r
+ fhPtHbpUeCharged->Fill(ptTrig,cosi);\r
+ if(DoEventSelect()){\r
+ for(Int_t im=0; im<GetMultiBin(); im++){\r
+ if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))\r
+ fhTrigUeCorr[im]->Fill(ptTrig,xE);\r
+ }\r
+ } //multiplicity events selection\r
+ \r
+ } //UE study\r
+ \r
+ if(fPi0Trigger){\r
+ if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
+ fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaphiDecay1);\r
+ fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaphiDecay2);\r
+ if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);\r
+ if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )\r
+ fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1); \r
+ if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )\r
+ fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);\r
+ if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);\r
+ } //index of decay photons found\r
+ } //make decay-hadron correlation \r
+ \r
+ //several UE calculation \r
+ if(fMakeSeveralUE){\r
+ if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){ \r
+ fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);\r
+ fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);\r
+ fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);\r
+ }\r
+ if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){ \r
+ fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);\r
+ fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);\r
+ fhPtHbpUeRightCharged->Fill(ptTrig,cosi);\r
+ \r
+ }\r
+ } //several UE calculation\r
+ \r
+ } //Fill histogram \r
+ else{\r
+ nrefs++;\r
+ if(nrefs==1){\r
+ reftracks = new TObjArray(0);\r
+ reftracks->SetName(GetAODObjArrayName()+"Tracks");\r
+ reftracks->SetOwner(kFALSE);\r
+ }\r
+ reftracks->Add(track);\r
+ }//aod particle loop\r
+ }// track loop\r
+ \r
+ //Fill AOD with reference tracks, if not filling histograms\r
+ if(!bFillHisto && reftracks) {\r
+ aodParticle->AddObjArray(reftracks);\r
+ }\r
+ \r
+ //delete reftracks;\r
+ \r
+} \r
+\r
+\r
+\r
+//____________________________________________________________________________\r
+//void AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector) \r
+//{ \r
+// // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed\r
+// if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");\r
+// \r
+// if(!NewOutputAOD()){\r
+// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");\r
+// abort();\r
+// }\r
+// \r
+// Double_t phiTrig = aodParticle->Phi();\r
+// Int_t tag = 0;\r
+// TLorentzVector gammai;\r
+// TLorentzVector gammaj;\r
+// \r
+// //Get vertex for photon momentum calculation\r
+// \r
+// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) \r
+// {\r
+// for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {\r
+// if (!GetMixedEvent()) \r
+// GetReader()->GetVertex(GetVertex(iev));\r
+// else \r
+// GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); \r
+// }\r
+// }\r
+// Double_t vertex2[] = {0.0,0.0,0.0} ; //vertex of second input aod\r
+// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) \r
+// {\r
+// if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);\r
+// }\r
+// \r
+// //Cluster loop, select pairs with good pt, phi and fill AODs or histograms\r
+// //Int_t iEvent= GetReader()->GetEventNumber() ;\r
+// Int_t nclus = pl->GetEntriesFast();\r
+// for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){\r
+// AliVCluster * calo = (AliVCluster *) (pl->At(iclus)) ;\r
+// \r
+// Int_t evtIndex1 = 0 ; \r
+// if (GetMixedEvent()) {\r
+// evtIndex1=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; \r
+// }\r
+// \r
+//\r
+// //Input from second AOD?\r
+// Int_t inputi = 0;\r
+// if (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= iclus) \r
+// inputi = 1 ;\r
+// else if(aodParticle->GetDetector() == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= iclus) \r
+// inputi = 1;\r
+// \r
+// //Cluster selection, not charged, with photon or pi0 id and in fiducial cut\r
+// //FIXME\r
+// Int_t pdg=0;\r
+// //if (inputi == 0 && !SelectCluster(calo, GetVertex(evtIndex1), gammai, pdg)) \r
+// continue ;\r
+// //MEFIX\r
+// else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg)) \r
+// continue ;\r
+// \r
+// if(GetDebug() > 2)\r
+// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral cluster in %s: pt %f, phi %f, phi trigger %f. Cuts: delta phi min %2.2f, max %2.2f, pT min %2.2f \n",\r
+// detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());\r
+// \r
+// //2 gamma overlapped, found with PID\r
+// if(pdg == AliCaloPID::kPi0){ \r
+// \r
+// //Select only hadrons in pt range\r
+// if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) \r
+// continue ;\r
+// \r
+// //Selection within angular range\r
+// Float_t phi = gammai.Phi();\r
+// if(phi < 0) phi+=TMath::TwoPi();\r
+// //Float_t deltaphi = TMath::Abs(phiTrig-phi);\r
+// //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;\r
+// \r
+// AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);\r
+// //pi0.SetLabel(calo->GetLabel());\r
+// pi0.SetPdg(AliCaloPID::kPi0);\r
+// pi0.SetDetector(detector);\r
+// \r
+// if(IsDataMC()){\r
+// pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(),inputi));\r
+// if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());\r
+// }//Work with stack also \r
+// //Set the indeces of the original caloclusters \r
+// pi0.SetCaloLabel(calo->GetID(),-1);\r
+// AddAODParticle(pi0);\r
+// \r
+// if(GetDebug() > 2) \r
+// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());\r
+// \r
+// }// pdg = 111\r
+// \r
+// //Make invariant mass analysis\r
+// else if(pdg == AliCaloPID::kPhoton){ \r
+// //Search the photon companion in case it comes from a Pi0 decay\r
+// //Apply several cuts to select the good pair;\r
+// for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){\r
+// AliVCluster * calo2 = (AliVCluster *) (pl->At(jclus)) ;\r
+// Int_t evtIndex2 = 0 ; \r
+// if (GetMixedEvent()) {\r
+// evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; \r
+// }\r
+// if (GetMixedEvent() && (evtIndex1 == evtIndex2))\r
+// continue ;\r
+// \r
+// //Input from second AOD?\r
+// Int_t inputj = 0;\r
+// if (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= jclus) \r
+// inputj = 1;\r
+// else if(aodParticle->GetDetector() == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= jclus) \r
+// inputj = 1;\r
+// \r
+// //Cluster selection, not charged with photon or pi0 id and in fiducial cut\r
+// Int_t pdgj=0;\r
+// //FIXME\r
+// //if (inputj == 0 && !SelectCluster(calo2, GetVertex(evtIndex2), gammaj, pdgj)) \r
+// continue ;\r
+// //MEFIX\r
+//\r
+// else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj)) \r
+// continue ;\r
+// //FIXME\r
+// //if(!SelectCluster(calo2,GetVertex(evtIndex2), gammaj, pdgj)) \r
+// //MEFIX\r
+// continue ;\r
+// \r
+// if(pdgj == AliCaloPID::kPhoton ){\r
+// \r
+// if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) \r
+// continue ;\r
+// \r
+// //Selection within angular range\r
+// Float_t phi = (gammai+gammaj).Phi();\r
+// if(phi < 0) phi+=TMath::TwoPi();\r
+// //Float_t deltaphi = TMath::Abs(phiTrig-phi);\r
+// //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;\r
+// \r
+// //Select good pair (aperture and invariant mass)\r
+// if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){\r
+// \r
+// if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",\r
+// (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());\r
+// \r
+// TLorentzVector pi0mom = gammai+gammaj;\r
+// AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);\r
+// //pi0.SetLabel(calo->GetLabel());\r
+// pi0.SetPdg(AliCaloPID::kPi0);\r
+// pi0.SetDetector(detector); \r
+// if(IsDataMC()){\r
+// //Check origin of the candidates\r
+// \r
+// Int_t label1 = calo->GetLabel();\r
+// Int_t label2 = calo2->GetLabel();\r
+// Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);\r
+// Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);\r
+// \r
+// if(GetDebug() > 0) \r
+// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
+// if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){\r
+// \r
+// //Check if pi0 mother is the same\r
+// if(GetReader()->ReadStack()){ \r
+// TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree\r
+// label1 = mother1->GetFirstMother();\r
+// //mother1 = GetMCStack()->Particle(label1);//pi0\r
+// \r
+// TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree\r
+// label2 = mother2->GetFirstMother();\r
+// //mother2 = GetMCStack()->Particle(label2);//pi0\r
+// }\r
+// else if(GetReader()->ReadAODMCParticles()){\r
+// AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree\r
+// label1 = mother1->GetMother();\r
+// //mother1 = GetMCStack()->Particle(label1);//pi0\r
+// AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree\r
+// label2 = mother2->GetMother();\r
+// //mother2 = GetMCStack()->Particle(label2);//pi0\r
+// }\r
+// \r
+// //printf("mother1 %d, mother2 %d\n",label1,label2);\r
+// if(label1 == label2)\r
+// GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);\r
+// }\r
+// }//Work with mc information also \r
+// pi0.SetTag(tag);\r
+// //Set the indeces of the original caloclusters \r
+// pi0.SetCaloLabel(calo->GetID(), calo2->GetID());\r
+// AddAODParticle(pi0);\r
+// \r
+// \r
+// }//Pair selected\r
+// }//if pair of gammas\r
+// }//2nd loop\r
+// }// if pdg = 22\r
+// }//1st loop\r
+// \r
+// if(GetDebug() > 1) \r
+// printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());\r
+//}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, TObjArray* pi0list, const Bool_t bFillHisto) \r
+{ \r
+ // Neutral Pion Correlation Analysis\r
+ if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());\r
+ \r
+ Int_t evtIndex11 = 0 ; \r
+ Int_t evtIndex12 = 0 ; \r
+ if (GetMixedEvent()) {\r
+ evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;\r
+ evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ; \r
+ } \r
+ \r
+ Double_t pt = -100.;\r
+ Double_t px = -100.;\r
+ Double_t py = -100.;\r
+ Double_t rat = -100.; \r
+ Double_t phi = -100.;\r
+ Double_t eta = -100.;\r
+ Double_t xE = -100.; \r
+ Double_t cosi = -100.; \r
+ \r
+ Double_t ptTrig = aodParticle->Pt();\r
+ Double_t phiTrig = aodParticle->Phi();\r
+ Double_t etaTrig = aodParticle->Eta();\r
+ Double_t pxTrig = aodParticle->Px();\r
+ Double_t pyTrig = aodParticle->Py();\r
+ \r
+ \r
+ Int_t indexPhoton1 = -1 ;\r
+ Int_t indexPhoton2 = -1 ; \r
+ Double_t ptDecay1 = 0. ;\r
+ Double_t pxDecay1 = 0. ;\r
+ Double_t pyDecay1 = 0. ;\r
+ Double_t phiDecay1 = 0. ;\r
+ Double_t ptDecay2 = 0. ;\r
+ Double_t pxDecay2 = 0. ;\r
+ Double_t pyDecay2 = 0. ;\r
+ Double_t phiDecay2 = 0. ;\r
+ \r
+ Double_t ratDecay1 = -100.; \r
+ Double_t ratDecay2 = -100.; \r
+ Float_t deltaphi = -100. ;\r
+ Float_t deltaphiDecay1 = -100. ;\r
+ Float_t deltaphiDecay2 = -100. ;\r
+ TObjArray * clusters = 0x0 ; \r
+ TLorentzVector photonMom ; \r
+ if(fPi0Trigger){\r
+ indexPhoton1 = aodParticle->GetCaloLabel (0);\r
+ indexPhoton2 = aodParticle->GetCaloLabel (1);\r
+ if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);\r
+ \r
+ if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
+ if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;\r
+ else clusters = GetPHOSClusters() ;\r
+ for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){\r
+ AliVCluster * photon = (AliVCluster*) (clusters->At(iclus)); \r
+ photon->GetMomentum(photonMom,GetVertex(0)) ;\r
+ if(photon->GetID()==indexPhoton1) {\r
+ ptDecay1 = photonMom.Pt();\r
+ pxDecay1 = photonMom.Px();\r
+ pyDecay1 = photonMom.Py();\r
+ phiDecay1 = photonMom.Phi();\r
+ }\r
+ if(photon->GetID()==indexPhoton2) {\r
+ ptDecay2 = photonMom.Pt();\r
+ pxDecay2 = photonMom.Px();\r
+ pyDecay2 = photonMom.Py();\r
+ phiDecay2 = photonMom.Phi();\r
+ } \r
+ if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);\r
+ } //photonAOD loop \r
+ } //index of decay photons found\r
+ if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);\r
+ } //make decay-hadron correlation\r
+ \r
+ TObjArray * refpi0 =0x0;\r
+ Int_t nrefs = 0;\r
+ \r
+ //Loop on stored AOD pi0\r
+ Int_t naod = pi0list->GetEntriesFast();\r
+ if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - aod branch entries %d\n", naod);\r
+ for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+ AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (pi0list->At(iaod));\r
+ \r
+ Int_t evtIndex2 = 0 ; \r
+ Int_t evtIndex3 = 0 ; \r
+ if (GetMixedEvent()) {\r
+ evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;\r
+ evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;\r
+ \r
+ if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events\r
+ continue ; \r
+ } \r
+ \r
+ //Int_t pdg = pi0->GetPdg();\r
+ //if(pdg != AliCaloPID::kPi0) continue; \r
+ \r
+ pt = pi0->Pt();\r
+ px = pi0->Px();\r
+ py = pi0->Py(); \r
+ if(pt < GetMinPt() || pt > GetMaxPt()) continue ;\r
+ //jumped out this event if near side associated partile pt larger than trigger\r
+ if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2()) break ;\r
+\r
+ //Selection within angular range\r
+ phi = pi0->Phi();\r
+ //Float_t deltaphi = TMath::Abs(phiTrig-phi);\r
+ //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;\r
+ \r
+ if(bFillHisto){\r
+ \r
+ deltaphi = phiTrig-phi;\r
+ if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();\r
+ if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();\r
+ \r
+ rat = pt/ptTrig ;\r
+ phi = pi0->Phi() ;\r
+ eta = pi0->Eta() ;\r
+ xE = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);\r
+ if(xE <0.)xE =-xE;\r
+ cosi = TMath::Log(1/xE);\r
+ \r
+ if(fPi0Trigger){\r
+ if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
+ if(ptDecay1) ratDecay1 = pt/ptDecay1 ;\r
+ if(ptDecay2) ratDecay2 = pt/ptDecay2 ; \r
+ deltaphiDecay1 = phiDecay1-phi;\r
+ deltaphiDecay2 = phiDecay2-phi;\r
+ if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();\r
+ if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();\r
+ if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();\r
+ if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi(); \r
+ fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaphiDecay1);\r
+ fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaphiDecay2);\r
+ if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);\r
+ if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )\r
+ fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1); \r
+ if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )\r
+ fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);\r
+ if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);\r
+ }\r
+ } //do decay-hadron correlation\r
+ \r
+ fhEtaNeutral->Fill(pt,eta);\r
+ fhPhiNeutral->Fill(pt,phi);\r
+ fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);\r
+ fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);\r
+ fhDeltaPhiDeltaEtaNeutral->Fill(deltaphi,etaTrig-eta);\r
+ \r
+ //delta phi cut for correlation\r
+ if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {\r
+ fhDeltaPhiNeutralPt->Fill(pt,deltaphi);\r
+ fhPtImbalanceNeutral->Fill(ptTrig,rat); \r
+ fhPtHbpNeutral->Fill(ptTrig,cosi); \r
+ }\r
+ else {\r
+ fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);\r
+ fhPtImbalanceUeNeutral->Fill(ptTrig,rat);\r
+ fhPtHbpUeNeutral->Fill(ptTrig,cosi); \r
+ }\r
+ //several UE calculation \r
+ if(fMakeSeveralUE){\r
+ if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){ \r
+ fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);\r
+ fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);\r
+ fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);\r
+ }\r
+ if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){ \r
+ fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);\r
+ fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);\r
+ fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);\r
+ }\r
+ } //several UE calculation\r
+ }\r
+ else{\r
+ nrefs++;\r
+ if(nrefs==1){\r
+ refpi0 = new TObjArray(0);\r
+ refpi0->SetName(GetAODObjArrayName()+"Pi0s");\r
+ refpi0->SetOwner(kFALSE);\r
+ }\r
+ refpi0->Add(pi0);\r
+ }//put references in trigger AOD \r
+ \r
+ //if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);\r
+ \r
+ }//loop\r
+}\r
+ \r
+\r
+//____________________________________________________________________________\r
+//Bool_t AliAnaParticleHadronCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {\r
+// //Select cluster depending on its pid and acceptance selections\r
+// \r
+// //Skip matched clusters with tracks\r
+// if(IsTrackMatched(calo)) return kFALSE;\r
+// \r
+// TString detector = "";\r
+// if (calo->IsPHOS()) detector= "PHOS";\r
+// else if(calo->IsEMCAL()) detector= "EMCAL";\r
+// \r
+// //Check PID\r
+// calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line\r
+// pdg = AliCaloPID::kPhoton; \r
+// if(IsCaloPIDOn()){\r
+// //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)\r
+// //or redo PID, recommended option for EMCal.\r
+// \r
+// if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )\r
+// pdg = GetCaloPID()->GetPdg(detector,calo->GetPID(),mom.E());//PID with weights\r
+// else\r
+// pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated\r
+// \r
+// if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);\r
+// \r
+// //If it does not pass pid, skip\r
+// if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {\r
+// return kFALSE ;\r
+// }\r
+// }//PID on\r
+// \r
+// //Check acceptance selection\r
+// if(IsFiducialCutOn()){\r
+// Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,detector) ;\r
+// if(! in ) return kFALSE ;\r
+// }\r
+// \r
+// if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);\r
+// \r
+// return kTRUE;\r
+// \r
+//}\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff