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