git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrected bug in cluster selection with PID
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleHadronCorrelation.cxx
CommitLineData
1c5acb87 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id: $ */
16
17//_________________________________________________________________________
18// Class for the analysis of particle - hadron correlations
19// Particle (for example direct gamma) must be found in a previous analysis
20//-- Author: Gustavo Conesa (LNF-INFN)
21//////////////////////////////////////////////////////////////////////////////
22
23
24// --- ROOT system ---
25#include "TH2F.h"
477d6cee 26#include "TClonesArray.h"
9415d854 27#include "TClass.h"
1c5acb87 28
29//---- ANALYSIS system ----
1c5acb87 30#include "AliNeutralMesonSelection.h"
31#include "AliAnaParticleHadronCorrelation.h"
32#include "AliCaloTrackReader.h"
33#include "AliCaloPID.h"
34#include "AliAODPWG4ParticleCorrelation.h"
477d6cee 35#include "AliFidutialCut.h"
36#include "AliAODTrack.h"
37#include "AliAODCaloCluster.h"
9415d854 38#include "AliMCAnalysisUtils.h"
39#include "TParticle.h"
40#include "AliStack.h"
1c5acb87 41
42ClassImp(AliAnaParticleHadronCorrelation)
43
44
45//____________________________________________________________________________
46 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation() :
47 AliAnaPartCorrBaseClass(),
6639984f 48 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
477d6cee 49 fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
1c5acb87 50 fhDeltaPhiCharged(0), fhDeltaPhiNeutral(0),
51 fhDeltaEtaCharged(0), fhDeltaEtaNeutral(0),
52 fhDeltaPhiChargedPt(0), fhDeltaPhiNeutralPt(0),
53 fhPtImbalanceNeutral(0), fhPtImbalanceCharged(0)
54{
55 //Default Ctor
477d6cee 56
1c5acb87 57 //Initialize parameters
58 InitParameters();
59}
60
61//____________________________________________________________________________
62AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(const AliAnaParticleHadronCorrelation & g) :
63 AliAnaPartCorrBaseClass(g),
64 fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
6639984f 65 fSelectIsolated(g.fSelectIsolated),
1c5acb87 66 fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral),
67 fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral),
68 fhDeltaPhiCharged(g.fhDeltaPhiCharged),
69 fhDeltaPhiNeutral(g.fhDeltaPhiNeutral),
70 fhDeltaEtaCharged(g.fhDeltaEtaCharged),
71 fhDeltaEtaNeutral(g.fhDeltaEtaNeutral),
72 fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt),
73 fhDeltaPhiNeutralPt(g.fhDeltaPhiNeutralPt),
74 fhPtImbalanceNeutral(g.fhPtImbalanceNeutral),
75 fhPtImbalanceCharged(g.fhPtImbalanceCharged)
76{
77 // cpy ctor
477d6cee 78
1c5acb87 79}
80
81//_________________________________________________________________________
82AliAnaParticleHadronCorrelation & AliAnaParticleHadronCorrelation::operator = (const AliAnaParticleHadronCorrelation & source)
83{
84 // assignment operator
477d6cee 85
1c5acb87 86 if(this == &source)return *this;
87 ((AliAnaPartCorrBaseClass *)this)->operator=(source);
88
89 fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
90 fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
6639984f 91 fSelectIsolated = source.fSelectIsolated ;
477d6cee 92
1c5acb87 93 fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ;
94 fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
95 fhDeltaPhiCharged = source.fhDeltaPhiCharged ;
96 fhDeltaPhiNeutral = source.fhDeltaPhiNeutral ;
97 fhDeltaPhiNeutralPt = source.fhDeltaPhiNeutralPt ;
98 fhDeltaEtaCharged = source.fhDeltaEtaCharged ;
99 fhDeltaEtaNeutral = source.fhDeltaEtaNeutral ;
100 fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
477d6cee 101
1c5acb87 102 fhPtImbalanceNeutral = source.fhPtImbalanceNeutral ;
103 fhPtImbalanceCharged = source.fhPtImbalanceCharged ;
104
105 return *this;
106
107}
108
109//________________________________________________________________________
110TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
111{
477d6cee 112
113 // Create histograms to be saved in output file and
114 // store them in fOutputContainer
115 TList * outputContainer = new TList() ;
116 outputContainer->SetName("CorrelationHistos") ;
117
118 Int_t nptbins = GetHistoNPtBins();
119 Int_t nphibins = GetHistoNPhiBins();
120 Int_t netabins = GetHistoNEtaBins();
121 Float_t ptmax = GetHistoPtMax();
122 Float_t phimax = GetHistoPhiMax();
123 Float_t etamax = GetHistoEtaMax();
124 Float_t ptmin = GetHistoPtMin();
125 Float_t phimin = GetHistoPhiMin();
126 Float_t etamin = GetHistoEtaMin();
127
128 //Correlation with charged hadrons
129 if(GetReader()->IsCTSSwitchedOn()) {
130 fhPhiCharged = new TH2F
131 ("PhiCharged","#phi_{h^{#pm}} vs p_{T trigger}",
132 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
133 fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
134 fhPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
135
136 fhEtaCharged = new TH2F
137 ("EtaCharged","#eta_{h^{#pm}} vs p_{T trigger}",
138 nptbins,ptmin,ptmax,netabins,etamin,etamax);
139 fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
140 fhEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
141
142 fhDeltaPhiCharged = new TH2F
143 ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
144 nptbins,ptmin,ptmax,200,0,6.4);
145 fhDeltaPhiCharged->SetYTitle("#Delta #phi");
146 fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
147
148 fhDeltaPhiChargedPt = new TH2F
9415d854 149 ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
477d6cee 150 nptbins,ptmin,ptmax,200,0,6.4);
151 fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
152 fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
153
154 fhDeltaEtaCharged = new TH2F
155 ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
156 nptbins,ptmin,ptmax,200,-2,2);
157 fhDeltaEtaCharged->SetYTitle("#Delta #eta");
158 fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
159
160 fhPtImbalanceCharged =
161 new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
162 nptbins,ptmin,ptmax,1000,0.,1.2);
163 fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
164 fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
165
166 outputContainer->Add(fhPhiCharged) ;
167 outputContainer->Add(fhEtaCharged) ;
168 outputContainer->Add(fhDeltaPhiCharged) ;
169 outputContainer->Add(fhDeltaEtaCharged) ;
170 outputContainer->Add(fhPtImbalanceCharged) ;
171 outputContainer->Add(fhDeltaPhiChargedPt) ;
172 } //Correlation with charged hadrons
173
174 //Correlation with neutral hadrons
175 if(GetReader()->IsEMCALSwitchedOn() || GetReader()->IsPHOSSwitchedOn()){
176
177 fhPhiNeutral = new TH2F
178 ("PhiNeutral","#phi_{#pi^{0}} vs p_{T trigger}",
179 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
180 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
181 fhPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
182
183 fhEtaNeutral = new TH2F
184 ("EtaNeutral","#eta_{#pi^{0}} vs p_{T trigger}",
185 nptbins,ptmin,ptmax,netabins,etamin,etamax);
186 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
187 fhEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
188
189 fhDeltaPhiNeutral = new TH2F
190 ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
191 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
192 fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
193 fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
194
195 fhDeltaPhiNeutralPt = new TH2F
196 ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
197 nptbins,ptmin,ptmax,200,0,6.4);
198 fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
199 fhDeltaPhiNeutralPt->SetXTitle("p_{T trigger} (GeV/c)");
200
201 fhDeltaEtaNeutral = new TH2F
202 ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
203 nptbins,ptmin,ptmax,200,-2,2);
204 fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
205 fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
206
207 fhPtImbalanceNeutral =
208 new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
209 nptbins,ptmin,ptmax,1000,0.,1.2);
210 fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
211 fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
212
213 outputContainer->Add(fhPhiNeutral) ;
214 outputContainer->Add(fhEtaNeutral) ;
215 outputContainer->Add(fhDeltaPhiNeutral) ;
216 outputContainer->Add(fhDeltaEtaNeutral) ;
217 outputContainer->Add(fhPtImbalanceNeutral) ;
218
219
220 //Keep neutral meson selection histograms if requiered
221 //Setting done in AliNeutralMesonSelection
222 if(GetNeutralMesonSelection()){
223 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
224 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
225 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
226 }
a3aebfff 227
477d6cee 228 }//Correlation with neutral hadrons
229
230 return outputContainer;
a3aebfff 231
1c5acb87 232}
233
477d6cee 234//____________________________________________________________________________
1c5acb87 235void AliAnaParticleHadronCorrelation::InitParameters()
236{
477d6cee 237
1c5acb87 238 //Initialize the parameters of the analysis.
a3aebfff 239 SetInputAODName("PWG4Particle");
240 SetAODRefArrayName("Hadrons");
241 AddToHistogramsName("AnaHadronCorr_");
242
9415d854 243 //Correlation with neutrals
244 //SetOutputAODClassName("AliAODPWG4Particle");
245 //SetOutputAODName("Pi0Correlated");
246
1c5acb87 247 SetPtCutRange(2,300);
248 fDeltaPhiMinCut = 1.5 ;
249 fDeltaPhiMaxCut = 4.5 ;
6639984f 250 fSelectIsolated = kFALSE;
1c5acb87 251}
252
253//__________________________________________________________________
254void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
255{
256
257 //Print some relevant parameters set for the analysis
258 if(! opt)
259 return;
260
a3aebfff 261 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
262 AliAnaPartCorrBaseClass::Print(" ");
263
1c5acb87 264 printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ;
265 printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;
6639984f 266 printf("Isolated Trigger? %d\n", fSelectIsolated) ;
1c5acb87 267}
268
269//____________________________________________________________________________
270void AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
271{
477d6cee 272 //Particle-Hadron Correlation Analysis, fill AODs
273
274 if(!GetInputAODBranch()){
a3aebfff 275 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, ABORT \n",GetInputAODName().Data());
477d6cee 276 abort();
277 }
9415d854 278
279 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
280 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());
281 abort();
282 }
283
477d6cee 284 if(GetDebug() > 1){
a3aebfff 285 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
286 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
287 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
288 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
289 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
477d6cee 290 }
291
292 //Loop on stored AOD particles, trigger
293 Int_t naod = GetInputAODBranch()->GetEntriesFast();
294 for(Int_t iaod = 0; iaod < naod ; iaod++){
295 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
9415d854 296
297 //Make correlation with charged hadrons
477d6cee 298 if(GetReader()->IsCTSSwitchedOn() )
a3aebfff 299 MakeChargedCorrelation(particle, GetAODCTS(),kFALSE);
477d6cee 300
301 //Make correlation with neutral pions
302 //Trigger particle in PHOS, correlation with EMCAL
303 if(particle->GetDetector()=="PHOS" && GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0)
9415d854 304 MakeNeutralCorrelationFillAOD(particle, GetAODEMCAL(),"EMCAL");
477d6cee 305 //Trigger particle in EMCAL, correlation with PHOS
306 else if(particle->GetDetector()=="EMCAL" && GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0)
9415d854 307 MakeNeutralCorrelationFillAOD(particle, GetAODPHOS(),"PHOS");
477d6cee 308 //Trigger particle in CTS, correlation with PHOS, EMCAL and CTS
309 else if(particle->GetDetector()=="CTS" ){
310 if(GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0)
9415d854 311 MakeNeutralCorrelationFillAOD(particle, GetAODPHOS(),"PHOS");
477d6cee 312 if(GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0)
9415d854 313 MakeNeutralCorrelationFillAOD(particle, GetAODEMCAL(),"EMCAL");
477d6cee 314 }
a3aebfff 315
316
477d6cee 317 }//Aod branch loop
318
a3aebfff 319 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
477d6cee 320
1c5acb87 321}
322
323//____________________________________________________________________________
324void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
325{
477d6cee 326 //Particle-Hadron Correlation Analysis, fill histograms
327
328 if(!GetInputAODBranch()){
a3aebfff 329 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, ABORT \n",GetInputAODName().Data());
477d6cee 330 abort();
331 }
9415d854 332
477d6cee 333 if(GetDebug() > 1){
a3aebfff 334 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
335 printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
477d6cee 336 }
9415d854 337
477d6cee 338 //Loop on stored AOD particles
339 Int_t naod = GetInputAODBranch()->GetEntriesFast();
9415d854 340 for(Int_t iaod = 0; iaod < naod ; iaod++){
341 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
342
343 //check if the particle is isolated or if we want to take the isolation into account
344 if(OnlyIsolated() && !particle->IsIsolated()) continue;
345
346 //Make correlation with charged hadrons
347 TRefArray * reftracks = particle->GetRefArray(GetAODRefArrayName()+"Tracks");
348 if(reftracks){
349 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Particle %d, In Track Refs entries %d\n", iaod, reftracks->GetEntriesFast());
350 if(reftracks->GetEntriesFast() > 0) MakeChargedCorrelation(particle, reftracks,kTRUE);
351 }
352
353 //Make correlation with neutral pions
354 if(GetOutputAODBranch() && GetOutputAODBranch()->GetEntriesFast() > 0){
355 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Particle %d, In Cluster Refs entries %d\n",iaod, GetOutputAODBranch()->GetEntriesFast());
356 MakeNeutralCorrelationFillHistograms(particle);
357 }
358
477d6cee 359 }//Aod branch loop
360
a3aebfff 361 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
477d6cee 362
1c5acb87 363}
364
365//____________________________________________________________________________
a3aebfff 366void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TRefArray* pl, const Bool_t bFillHisto)
1c5acb87 367{
368 // Charged Hadron Correlation Analysis
a3aebfff 369 if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
1c5acb87 370
371 Double_t ptTrig = aodParticle->Pt();
372 Double_t phiTrig = aodParticle->Phi();
373 Double_t pt = -100.;
374 Double_t rat = -100.;
375 Double_t phi = -100. ;
376 Double_t eta = -100. ;
377 Double_t p[3];
477d6cee 378 Bool_t first=kTRUE;
1c5acb87 379
a3aebfff 380 TRefArray * reftracks =0x0;
381 if(!bFillHisto)
382 reftracks = new TRefArray;
383
384
1c5acb87 385 //Track loop, select tracks with good pt, phi and fill AODs or histograms
386 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
387 AliAODTrack * track = (AliAODTrack *) (pl->At(ipr)) ;
388 track->GetPxPyPz(p) ;
389 TLorentzVector mom(p[0],p[1],p[2],0);
9415d854 390 pt = mom.Pt();
1c5acb87 391 eta = mom.Eta();
392 phi = mom.Phi() ;
9415d854 393 if(phi < 0) phi+=TMath::TwoPi();
1c5acb87 394 rat = pt/ptTrig ;
477d6cee 395
1c5acb87 396 if(IsFidutialCutOn()){
397 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
398 if(! in ) continue ;
399 }
400
401 //Select only hadrons in pt range
402 if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
403
9415d854 404 //Selection within angular range
405 Float_t deltaphi = TMath::Abs(phiTrig-phi);
406 if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
407
1c5acb87 408 if(GetDebug() > 2)
9415d854 409 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",
410 pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());
411
1c5acb87 412 if(bFillHisto){
413 // Fill Histograms
414 fhEtaCharged->Fill(ptTrig,eta);
415 fhPhiCharged->Fill(ptTrig,phi);
416 fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
9415d854 417 fhDeltaPhiCharged->Fill(ptTrig, deltaphi);
418 fhDeltaPhiChargedPt->Fill(pt,deltaphi);
419 if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
420 fhPtImbalanceCharged->Fill(ptTrig,rat);
1c5acb87 421 }
422 else{
423 //Fill AODs
477d6cee 424
425 if(first) {
a3aebfff 426 new (reftracks) TRefArray(TProcessID::GetProcessWithUID(track));
477d6cee 427 first = kFALSE;
428 }
9415d854 429
a3aebfff 430 reftracks->Add(track);
1c5acb87 431
432 }//aod particle loop
433 }// track loop
9415d854 434
a3aebfff 435 //Fill AOD with reference tracks, if not filling histograms
436 if(!bFillHisto && reftracks->GetEntriesFast() > 0) {
437 reftracks->SetName(GetAODRefArrayName()+"Tracks");
438 aodParticle->AddRefArray(reftracks);
439 }
9415d854 440
1c5acb87 441}
9415d854 442
1c5acb87 443//____________________________________________________________________________
9415d854 444void AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation * aodParticle,TRefArray* pl, TString detector)
1c5acb87 445{
9415d854 446 // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed
447 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");
448
449 if(!NewOutputAOD()){
450 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, ABORT! \n");
451 abort();
452 }
1c5acb87 453
1c5acb87 454 Double_t phiTrig = aodParticle->Phi();
9415d854 455 Int_t tag = -1;
1c5acb87 456 TLorentzVector gammai;
457 TLorentzVector gammaj;
458
459 Double_t vertex[] = {0,0,0};
1c5acb87 460 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
461
462 //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
9415d854 463 //Int_t iEvent= GetReader()->GetEventNumber() ;
464 Int_t nclus = pl->GetEntriesFast();
465 for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){
1c5acb87 466 AliAODCaloCluster * calo = (AliAODCaloCluster *) (pl->At(iclus)) ;
9415d854 467
468 //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
469 Int_t pdg=0;
470 if(!SelectCluster(calo, vertex, gammai, pdg)) continue ;
12524a23 471
472 if(GetDebug() > 2)
9415d854 473 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",
474 detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
475
476 //2 gamma overlapped, found with PID
477 if(pdg == AliCaloPID::kPi0){
477d6cee 478
9415d854 479 //Select only hadrons in pt range
480 if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) continue ;
1c5acb87 481
9415d854 482 //Selection within angular range
483 Float_t phi = gammai.Phi();
484 if(phi < 0) phi+=TMath::TwoPi();
485 Float_t deltaphi = TMath::Abs(phiTrig-phi);
486 if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
477d6cee 487
9415d854 488 AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);
489 //pi0.SetLabel(calo->GetLabel(0));
490 pi0.SetPdg(AliCaloPID::kPi0);
491 pi0.SetDetector(detector);
492
493 if(IsDataMC()){
494 pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
495 if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());
496 }//Work with stack also
497 //Set the indeces of the original caloclusters
498 pi0.SetCaloLabel(calo->GetID(),-1);
499 AddAODParticle(pi0);
500
501 if(GetDebug() > 2)
502 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());
503
504 }// pdg = 111
505
506 //Make invariant mass analysis
507 else if(pdg == AliCaloPID::kPhoton){
508 //Search the photon companion in case it comes from a Pi0 decay
509 //Apply several cuts to select the good pair;
510 for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
511 AliAODCaloCluster * calo2 = (AliAODCaloCluster *) (pl->At(jclus)) ;
477d6cee 512
9415d854 513 //Cluster selection, not charged with photon or pi0 id and in fidutial cut
514 Int_t pdgj=0;
515 if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
477d6cee 516
9415d854 517 if(pdgj == AliCaloPID::kPhoton ){
518
519 if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) continue ;
1c5acb87 520
9415d854 521 //Selection within angular range
522 Float_t phi = (gammai+gammaj).Phi();
523 if(phi < 0) phi+=TMath::TwoPi();
524 Float_t deltaphi = TMath::Abs(phiTrig-phi);
525 if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
477d6cee 526
9415d854 527 //Select good pair (aperture and invariant mass)
528 if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
1c5acb87 529
9415d854 530 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",
531 (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
1c5acb87 532
9415d854 533 TLorentzVector pi0mom = gammai+gammaj;
534 AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);
535 //pi0.SetLabel(calo->GetLabel(0));
536 pi0.SetPdg(AliCaloPID::kPi0);
537 pi0.SetDetector(detector);
538 if(IsDataMC()){
539 //Check origin of the candidates
540 Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0), GetMCStack());
541 Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabel(0), GetMCStack());
477d6cee 542
9415d854 543 if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
544 if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){
545
546 //Check if pi0 mother is the same
547 Int_t label1 = calo->GetLabel(0);
548 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
549 label1 = mother1->GetFirstMother();
550 //mother1 = GetMCStack()->Particle(label1);//pi0
551
552 Int_t label2 = calo2->GetLabel(0);
553 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
554 label2 = mother2->GetFirstMother();
555 //mother2 = GetMCStack()->Particle(label2);//pi0
556
557 //printf("mother1 %d, mother2 %d\n",label1,label2);
558 if(label1 == label2)
559 tag = AliMCAnalysisUtils::kMCPi0;
477d6cee 560 }
9415d854 561 }//Work with stack also
562 pi0.SetTag(tag);
563 //Set the indeces of the original caloclusters
564 pi0.SetCaloLabel(calo->GetID(), calo2->GetID());
565 AddAODParticle(pi0);
566
567
568 }//Pair selected
569 }//if pair of gammas
570 }//2nd loop
571 }// if pdg = 22
1c5acb87 572 }//1st loop
573
9415d854 574 if(GetDebug() > 1)
575 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());
576}
577
578//____________________________________________________________________________
579void AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms(AliAODPWG4ParticleCorrelation * aodParticle)
580{
581 // Neutral Pion Correlation Analysis
582 if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistogramS() - Make trigger particle - pi0 correlation, %d pi0's \n",GetOutputAODBranch()->GetEntriesFast());
583
584 Double_t pt = -100.;
585 Double_t rat = -100.;
586 Double_t phi = -100.;
587 Double_t eta = -100.;
588 Double_t ptTrig = aodParticle->Pt();
589 Double_t phiTrig = aodParticle->Phi();
590 Double_t etaTrig = aodParticle->Eta();
591
592 if(!GetOutputAODBranch()){
593 printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
594 abort();
a3aebfff 595 }
9415d854 596
597 //Loop on stored AOD pi0
598 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
599 if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - aod branch entries %d\n", naod);
600 for(Int_t iaod = 0; iaod < naod ; iaod++){
601 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
602 Int_t pdg = pi0->GetPdg();
603
604 if(pdg != AliCaloPID::kPi0) continue;
605 pt = pi0->Pt();
606
607 if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
608
609 //Selection within angular range
610 phi = pi0->Phi();
611 Float_t deltaphi = TMath::Abs(phiTrig-phi);
612 if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
613
614 rat = pt/ptTrig ;
615 phi = pi0->Phi() ;
616 eta = pi0->Eta() ;
617
618 fhEtaNeutral->Fill(ptTrig,eta);
619 fhPhiNeutral->Fill(ptTrig,phi);
620 fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
621 fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);
622 fhDeltaPhiNeutralPt->Fill(pt,deltaphi);
623 fhPtImbalanceNeutral->Fill(ptTrig,rat);
624
625 if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
626
627 }//loop
1c5acb87 628}
629
9415d854 630
1c5acb87 631//____________________________________________________________________________
632Bool_t AliAnaParticleHadronCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) const {
9415d854 633 //Select cluster depending on its pid and acceptance selections
634
635 //Skip matched clusters with tracks
12524a23 636 if(calo->GetNTracksMatched() > 0) return kFALSE;
9415d854 637
12524a23 638 TString detector = "";
639 if(calo->IsPHOSCluster()) detector= "PHOS";
640 else if(calo->IsEMCALCluster()) detector= "EMCAL";
641
9415d854 642 //Check PID
643 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
644 pdg = AliCaloPID::kPhoton;
645 if(IsCaloPIDOn()){
646 //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
647 //or redo PID, recommended option for EMCal.
9415d854 648
649 if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
650 pdg = GetCaloPID()->GetPdg(detector,calo->PID(),mom.E());//PID with weights
651 else
652 pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
653
654 if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
655
656 //If it does not pass pid, skip
12524a23 657 if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {
9415d854 658 return kFALSE ;
659 }
12524a23 660 }//PID on
9415d854 661
662 //Check acceptance selection
663 if(IsFidutialCutOn()){
12524a23 664 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,detector) ;
665 if(! in ) return kFALSE ;
9415d854 666 }
667
668 if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
669
670 return kTRUE;
671
672}