]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaParticleIsolation.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Data are not reset when method WriteVertices is called. Added debug histograms
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleIsolation.cxx
CommitLineData
1a31a9ab 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 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: AliAnaParticleIsolation.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18// Class for analysis of particle isolation
19// Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)
20//
21// Class created from old AliPHOSGammaJet
22// (see AliRoot versions previous Release 4-09)
23//
24// -- Author: Gustavo Conesa (LNF-INFN)
25
26//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
27//////////////////////////////////////////////////////////////////////////////
28
29
30// --- ROOT system ---
31#include <TClonesArray.h>
32#include <TList.h>
33#include <TObjString.h>
34#include <TH2F.h>
35//#include <Riostream.h>
36#include <TClass.h>
37
38// --- Analysis system ---
39#include "AliAnaParticleIsolation.h"
40#include "AliCaloTrackReader.h"
41#include "AliIsolationCut.h"
42#include "AliNeutralMesonSelection.h"
43#include "AliAODPWG4ParticleCorrelation.h"
44#include "AliMCAnalysisUtils.h"
45#include "AliVTrack.h"
46#include "AliVCluster.h"
47
48ClassImp(AliAnaParticleIsolation)
49
50//____________________________________________________________________________
51 AliAnaParticleIsolation::AliAnaParticleIsolation() :
52 AliAnaPartCorrBaseClass(), fCalorimeter(""),
53 fReMakeIC(0), fMakeSeveralIC(0), fMakeInvMass(0),
54 fhPtIso(0),fhPhiIso(0),fhEtaIso(0), fhPtNoIso(0),fhPtInvMassDecayIso(0),fhPtInvMassDecayNoIso(0), fhConeSumPt(0), fhPtInCone(0),
55 fhFRConeSumPt(0), fhPtInFRCone(0),
56 //Several IC
57 fNCones(0),fNPtThresFrac(0), fConeSizes(), fPtThresholds(), fPtFractions(),
58 //MC
59 fhPtIsoPrompt(0),fhPhiIsoPrompt(0),fhEtaIsoPrompt(0),
60 fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(), fhPtSumIsolatedPrompt(),
61 fhPtIsoFragmentation(0),fhPhiIsoFragmentation(0),fhEtaIsoFragmentation(0),
62 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
63 fhPtIsoPi0Decay(0),fhPhiIsoPi0Decay(0),fhEtaIsoPi0Decay(0),
64 fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(), fhPtSumIsolatedPi0Decay(),
65 fhPtIsoOtherDecay(0),fhPhiIsoOtherDecay(0),fhEtaIsoOtherDecay(0),
66 fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(), fhPtSumIsolatedOtherDecay(),
67 fhPtIsoConversion(0),fhPhiIsoConversion(0),fhEtaIsoConversion(0),
68 fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(), fhPtSumIsolatedConversion(),
69 fhPtIsoUnknown(0),fhPhiIsoUnknown(0),fhEtaIsoUnknown(0),
70 fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(), fhPtSumIsolatedUnknown(),
71 fhPtNoIsoPi0Decay(0),fhPtNoIsoPrompt(0),fhPtIsoMCPhoton(0),fhPtNoIsoMCPhoton(0),
72 //Histograms settings
73 fHistoNPtSumBins(0), fHistoPtSumMax(0.), fHistoPtSumMin(0.),
74 fHistoNPtInConeBins(0), fHistoPtInConeMax(0.), fHistoPtInConeMin(0.)
75{
76 //default ctor
77
78 //Initialize parameters
79 InitParameters();
80
81 for(Int_t i = 0; i < 5 ; i++){
82 fConeSizes[i] = 0 ;
83 fhPtSumIsolated[i] = 0 ;
84
85 fhPtSumIsolatedPrompt[i] = 0 ;
86 fhPtSumIsolatedFragmentation[i] = 0 ;
87 fhPtSumIsolatedPi0Decay[i] = 0 ;
88 fhPtSumIsolatedOtherDecay[i] = 0 ;
89 fhPtSumIsolatedConversion[i] = 0 ;
90 fhPtSumIsolatedUnknown[i] = 0 ;
91
92 for(Int_t j = 0; j < 5 ; j++){
93 fhPtThresIsolated[i][j] = 0 ;
94 fhPtFracIsolated[i][j] = 0 ;
95
96 fhPtThresIsolatedPrompt[i][j] = 0 ;
97 fhPtThresIsolatedFragmentation[i][j] = 0 ;
98 fhPtThresIsolatedPi0Decay[i][j] = 0 ;
99 fhPtThresIsolatedOtherDecay[i][j] = 0 ;
100 fhPtThresIsolatedConversion[i][j] = 0 ;
101 fhPtThresIsolatedUnknown[i][j] = 0 ;
102
103 fhPtFracIsolatedPrompt[i][j] = 0 ;
104 fhPtFracIsolatedFragmentation[i][j] = 0 ;
105 fhPtFracIsolatedPi0Decay[i][j] = 0 ;
106 fhPtFracIsolatedOtherDecay[i][j] = 0 ;
107 fhPtFracIsolatedConversion[i][j] = 0 ;
108 fhPtFracIsolatedUnknown[i][j] = 0 ;
109
110 }
111 }
112
113 for(Int_t i = 0; i < 5 ; i++){
114 fPtFractions[i]= 0 ;
115 fPtThresholds[i]= 0 ;
116 }
117
118
119}
120
121//____________________________________________________________________________
122AliAnaParticleIsolation::~AliAnaParticleIsolation()
123{
124 //dtor
125 //do not delete histograms
126
127 //delete [] fConeSizes ;
128 //delete [] fPtThresholds ;
129 //delete [] fPtFractions ;
130
131}
132
133//_________________________________________________________________________
134Bool_t AliAnaParticleIsolation::CheckInvMass(const Int_t iaod, const AliAODPWG4Particle * part1)
135{
136 // Search if there is a companion decay photon to the candidate
137 // and discard it in such case
138 // Use it only if isolation candidates are photons
139 // Make sure that no selection on photon pt is done in the input aod photon list.
140 TLorentzVector mom1 = *(part1->Momentum());
141 TLorentzVector mom2 ;
142 for(Int_t jaod = 0; jaod < GetInputAODBranch()->GetEntriesFast(); jaod++){
143 if(iaod == jaod) continue ;
144 AliAODPWG4ParticleCorrelation * part2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jaod));
145 mom2 = *(part2->Momentum());
3bfcb597 146
1a31a9ab 147 //Select good pair (good phi, pt cuts, aperture and invariant mass)
3bfcb597 148
149 TString detector="";
150 if(part1->GetDetector()=="EMCAL" || part2->GetDetector()=="EMCAL" )
151 detector = "EMCAL";
152
153 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,detector)){
1a31a9ab 154 if(GetDebug() > 1)printf("AliAnaParticleIsolation::CheckInvMass() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
155 return kTRUE ;
156 }
157 }//loop
158
159 return kFALSE;
160}
161
162//________________________________________________________________________
163TObjString * AliAnaParticleIsolation::GetAnalysisCuts()
164{
165 //Save parameters used for analysis
166 TString parList ; //this will be list of parameters used for this analysis.
167 const Int_t buffersize = 255;
168 char onePar[buffersize] ;
169
170 snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
171 parList+=onePar ;
172 snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
173 parList+=onePar ;
174 snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
175 parList+=onePar ;
176 snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
177 parList+=onePar ;
178 snprintf(onePar, buffersize,"fMakeInvMass=%d (Flag for rejection of candidates with a pi0 inv mass pair) \n",fMakeInvMass) ;
179 parList+=onePar ;
180
181 if(fMakeSeveralIC){
182 snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
183 parList+=onePar ;
184 snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
185 parList+=onePar ;
186
187 for(Int_t icone = 0; icone < fNCones ; icone++){
188 snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
189 parList+=onePar ;
190 }
191 for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
192 snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
193 parList+=onePar ;
194 }
195 for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++){
196 snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
197 parList+=onePar ;
198 }
199 }
200
201 //Get parameters set in base class.
202 parList += GetBaseParametersList() ;
203
204 //Get parameters set in IC class.
205 if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
206
207 return new TObjString(parList) ;
208
209}
210
211//________________________________________________________________________
212TList * AliAnaParticleIsolation::GetCreateOutputObjects()
213{
214 // Create histograms to be saved in output file and
215 // store them in outputContainer
216 TList * outputContainer = new TList() ;
217 outputContainer->SetName("IsolatedParticleHistos") ;
218
219 Int_t nptbins = GetHistoPtBins();
220 Int_t nphibins = GetHistoPhiBins();
221 Int_t netabins = GetHistoEtaBins();
222 Float_t ptmax = GetHistoPtMax();
223 Float_t phimax = GetHistoPhiMax();
224 Float_t etamax = GetHistoEtaMax();
225 Float_t ptmin = GetHistoPtMin();
226 Float_t phimin = GetHistoPhiMin();
227 Float_t etamin = GetHistoEtaMin();
228
229 Int_t nptsumbins = fHistoNPtSumBins;
230 Float_t ptsummax = fHistoPtSumMax;
231 Float_t ptsummin = fHistoPtSumMin;
232 Int_t nptinconebins = fHistoNPtInConeBins;
233 Float_t ptinconemax = fHistoPtInConeMax;
234 Float_t ptinconemin = fHistoPtInConeMin;
235
236 if(!fMakeSeveralIC){
237
238 fhConeSumPt = new TH2F
239 ("hConePtSum","#Sigma p_{T} in isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
240 fhConeSumPt->SetYTitle("#Sigma p_{T}");
241 fhConeSumPt->SetXTitle("p_{T} (GeV/c)");
242 outputContainer->Add(fhConeSumPt) ;
243
244 fhPtInCone = new TH2F
245 ("hPtInCone","p_{T} in isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
246 fhPtInCone->SetYTitle("p_{T in cone} (GeV/c)");
247 fhPtInCone->SetXTitle("p_{T} (GeV/c)");
248 outputContainer->Add(fhPtInCone) ;
249
250 fhFRConeSumPt = new TH2F
251 ("hFRConePtSum","#Sigma p_{T} in the froward region isolation cone ",nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
252 fhFRConeSumPt->SetYTitle("#Sigma p_{T}");
253 fhFRConeSumPt->SetXTitle("p_{T} (GeV/c)");
254 outputContainer->Add(fhFRConeSumPt) ;
255
256 fhPtInFRCone = new TH2F
257 ("hPtInFRCone","p_{T} in froward region isolation cone ",nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
258 fhPtInFRCone->SetYTitle("p_{T in cone} (GeV/c)");
259 fhPtInFRCone->SetXTitle("p_{T} (GeV/c)");
260 outputContainer->Add(fhPtInFRCone) ;
261
732895a6 262 fhPtIso = new TH1F("hPt","Number of isolated particles",nptbins,ptmin,ptmax);
1a31a9ab 263 fhPtIso->SetYTitle("N");
264 fhPtIso->SetXTitle("p_{T}(GeV/c)");
265 outputContainer->Add(fhPtIso) ;
266
267 fhPhiIso = new TH2F
732895a6 268 ("hPhi","Number of isolated particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 269 fhPhiIso->SetYTitle("#phi");
270 fhPhiIso->SetXTitle("p_{T} (GeV/c)");
271 outputContainer->Add(fhPhiIso) ;
272
273 fhEtaIso = new TH2F
732895a6 274 ("hEta","Number of isolated particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 275 fhEtaIso->SetYTitle("#eta");
276 fhEtaIso->SetXTitle("p_{T} (GeV/c)");
277 outputContainer->Add(fhEtaIso) ;
278
279 fhPtNoIso = new TH1F("hPtNoIso","Number of not isolated leading particles",nptbins,ptmin,ptmax);
280 fhPtNoIso->SetYTitle("N");
281 fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
282 outputContainer->Add(fhPtNoIso) ;
283
732895a6 284 fhPtInvMassDecayIso = new TH1F("hPtInvMassDecayIso","Number of isolated #pi^{0} decay particles",nptbins,ptmin,ptmax);
1a31a9ab 285 fhPtNoIso->SetYTitle("N");
286 fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
287 outputContainer->Add(fhPtInvMassDecayIso) ;
288
289 fhPtInvMassDecayNoIso = new TH1F("hPtInvMassDecayNoIso","Number of not isolated leading pi0 decay particles",nptbins,ptmin,ptmax);
290 fhPtNoIso->SetYTitle("N");
291 fhPtNoIso->SetXTitle("p_{T}(GeV/c)");
292 outputContainer->Add(fhPtInvMassDecayNoIso) ;
293
294 if(IsDataMC()){
295
732895a6 296 fhPtIsoPrompt = new TH1F("hPtMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax);
1a31a9ab 297 fhPtIsoPrompt->SetYTitle("N");
298 fhPtIsoPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
299 outputContainer->Add(fhPtIsoPrompt) ;
300
301 fhPhiIsoPrompt = new TH2F
732895a6 302 ("hPhiMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 303 fhPhiIsoPrompt->SetYTitle("#phi");
304 fhPhiIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
305 outputContainer->Add(fhPhiIsoPrompt) ;
306
307 fhEtaIsoPrompt = new TH2F
732895a6 308 ("hEtaMCPrompt","Number of isolated prompt #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 309 fhEtaIsoPrompt->SetYTitle("#eta");
310 fhEtaIsoPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
311 outputContainer->Add(fhEtaIsoPrompt) ;
312
732895a6 313 fhPtIsoFragmentation = new TH1F("hPtMCFragmentation","Number of isolated #gamma",nptbins,ptmin,ptmax);
1a31a9ab 314 fhPtIsoFragmentation->SetYTitle("N");
315 fhPtIsoFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
316 outputContainer->Add(fhPtIsoFragmentation) ;
317
318 fhPhiIsoFragmentation = new TH2F
732895a6 319 ("hPhiMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 320 fhPhiIsoFragmentation->SetYTitle("#phi");
321 fhPhiIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
322 outputContainer->Add(fhPhiIsoFragmentation) ;
323
324 fhEtaIsoFragmentation = new TH2F
732895a6 325 ("hEtaMCFragmentation","Number of isolated fragmentation #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 326 fhEtaIsoFragmentation->SetYTitle("#eta");
327 fhEtaIsoFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
328 outputContainer->Add(fhEtaIsoFragmentation) ;
329
732895a6 330 fhPtIsoPi0Decay = new TH1F("hPtMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
1a31a9ab 331 fhPtIsoPi0Decay->SetYTitle("N");
332 fhPtIsoPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
333 outputContainer->Add(fhPtIsoPi0Decay) ;
334
335 fhPhiIsoPi0Decay = new TH2F
732895a6 336 ("hPhiMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 337 fhPhiIsoPi0Decay->SetYTitle("#phi");
338 fhPhiIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
339 outputContainer->Add(fhPhiIsoPi0Decay) ;
340
341 fhEtaIsoPi0Decay = new TH2F
732895a6 342 ("hEtaMCPi0Decay","Number of isolated #gamma from #pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 343 fhEtaIsoPi0Decay->SetYTitle("#eta");
344 fhEtaIsoPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
345 outputContainer->Add(fhEtaIsoPi0Decay) ;
346
732895a6 347 fhPtIsoOtherDecay = new TH1F("hPtMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax);
1a31a9ab 348 fhPtIsoOtherDecay->SetYTitle("N");
349 fhPtIsoOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
350 outputContainer->Add(fhPtIsoOtherDecay) ;
351
352 fhPhiIsoOtherDecay = new TH2F
732895a6 353 ("hPhiMCOtherDecay","Number of isolated #gamma from non-#pi^{0} decay",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 354 fhPhiIsoOtherDecay->SetYTitle("#phi");
355 fhPhiIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
356 outputContainer->Add(fhPhiIsoOtherDecay) ;
357
358 fhEtaIsoOtherDecay = new TH2F
732895a6 359 ("hEtaMCOtherDecay","Number of isolated #gamma non-#pi^{0} decay",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 360 fhEtaIsoOtherDecay->SetYTitle("#eta");
361 fhEtaIsoOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
362 outputContainer->Add(fhEtaIsoOtherDecay) ;
363
732895a6 364 fhPtIsoConversion = new TH1F("hPtMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax);
1a31a9ab 365 fhPtIsoConversion->SetYTitle("N");
366 fhPtIsoConversion->SetXTitle("p_{T #gamma}(GeV/c)");
367 outputContainer->Add(fhPtIsoConversion) ;
368
369 fhPhiIsoConversion = new TH2F
732895a6 370 ("hPhiMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 371 fhPhiIsoConversion->SetYTitle("#phi");
372 fhPhiIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
373 outputContainer->Add(fhPhiIsoConversion) ;
374
375 fhEtaIsoConversion = new TH2F
732895a6 376 ("hEtaMCConversion","Number of isolated converted #gamma",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 377 fhEtaIsoConversion->SetYTitle("#eta");
378 fhEtaIsoConversion->SetXTitle("p_{T #gamma} (GeV/c)");
379 outputContainer->Add(fhEtaIsoConversion) ;
380
732895a6 381 fhPtIsoUnknown = new TH1F("hPtMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax);
1a31a9ab 382 fhPtIsoUnknown->SetYTitle("N");
383 fhPtIsoUnknown->SetXTitle("p_{T}(GeV/c)");
384 outputContainer->Add(fhPtIsoUnknown) ;
385
386 fhPhiIsoUnknown = new TH2F
732895a6 387 ("hPhiMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1a31a9ab 388 fhPhiIsoUnknown->SetYTitle("#phi");
389 fhPhiIsoUnknown->SetXTitle("p_{T} (GeV/c)");
390 outputContainer->Add(fhPhiIsoUnknown) ;
391
392 fhEtaIsoUnknown = new TH2F
732895a6 393 ("hEtaMCUnknown","Number of isolated non-#gamma particles",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1a31a9ab 394 fhEtaIsoUnknown->SetYTitle("#eta");
395 fhEtaIsoUnknown->SetXTitle("p_{T} (GeV/c)");
396 outputContainer->Add(fhEtaIsoUnknown) ;
397
398 fhPtNoIsoPi0Decay = new TH1F
732895a6 399 ("hPtNoIsoPi0Decay","Number of not isolated leading #gamma from #pi^{0} decay",nptbins,ptmin,ptmax);
1a31a9ab 400 fhPtNoIsoPi0Decay->SetYTitle("N");
401 fhPtNoIsoPi0Decay->SetXTitle("p_{T} (GeV/c)");
402 outputContainer->Add(fhPtNoIsoPi0Decay) ;
403
404 fhPtNoIsoPrompt = new TH1F
405 ("hPtNoIsoPrompt","Number of not isolated leading prompt #gamma",nptbins,ptmin,ptmax);
406 fhPtNoIsoPrompt->SetYTitle("N");
407 fhPtNoIsoPrompt->SetXTitle("p_{T} (GeV/c)");
408 outputContainer->Add(fhPtNoIsoPrompt) ;
409
410 fhPtIsoMCPhoton = new TH1F
411 ("hPtIsoMCPhoton","Number of isolated leading #gamma",nptbins,ptmin,ptmax);
412 fhPtIsoMCPhoton->SetYTitle("N");
413 fhPtIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
414 outputContainer->Add(fhPtIsoMCPhoton) ;
415
416 fhPtNoIsoMCPhoton = new TH1F
732895a6 417 ("hPtNoIsoMCPhoton","Number of not isolated leading #gamma",nptbins,ptmin,ptmax);
1a31a9ab 418 fhPtNoIsoMCPhoton->SetYTitle("N");
419 fhPtNoIsoMCPhoton->SetXTitle("p_{T} (GeV/c)");
420 outputContainer->Add(fhPtNoIsoMCPhoton) ;
421 }//Histos with MC
422
423 }
424
425 if(fMakeSeveralIC){
426 const Int_t buffersize = 255;
427 char name[buffersize];
428 char title[buffersize];
429 for(Int_t icone = 0; icone<fNCones; icone++){
430 snprintf(name, buffersize,"hPtSum_Cone_%d",icone);
431 snprintf(title, buffersize,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
432 fhPtSumIsolated[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
433 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
434 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
435 outputContainer->Add(fhPtSumIsolated[icone]) ;
436
437 if(IsDataMC()){
438 snprintf(name, buffersize,"hPtSumPrompt_Cone_%d",icone);
439 snprintf(title, buffersize,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
440 fhPtSumIsolatedPrompt[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
441 fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
442 fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
443 outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
444
445 snprintf(name, buffersize,"hPtSumFragmentation_Cone_%d",icone);
446 snprintf(title, buffersize,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
447 fhPtSumIsolatedFragmentation[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
448 fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
449 fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
450 outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
451
452 snprintf(name, buffersize,"hPtSumPi0Decay_Cone_%d",icone);
453 snprintf(title, buffersize,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
454 fhPtSumIsolatedPi0Decay[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
455 fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
456 fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
457 outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
458
459 snprintf(name, buffersize,"hPtSumOtherDecay_Cone_%d",icone);
460 snprintf(title, buffersize,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
461 fhPtSumIsolatedOtherDecay[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
462 fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
463 fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
464 outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
465
466 snprintf(name, buffersize,"hPtSumConversion_Cone_%d",icone);
467 snprintf(title, buffersize,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
468 fhPtSumIsolatedConversion[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
469 fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
470 fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
471 outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
472
473 snprintf(name, buffersize,"hPtSumUnknown_Cone_%d",icone);
474 snprintf(title, buffersize,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
475 fhPtSumIsolatedUnknown[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
476 fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
477 fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
478 outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
479
480 }//Histos with MC
481
482 for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){
483 snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
484 snprintf(title, buffersize,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
485 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
486 fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
487 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
488
489 snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
490 snprintf(title, buffersize,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
491 fhPtFracIsolated[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
492 fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
493 outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
494
495 if(IsDataMC()){
496 snprintf(name, buffersize,"hPtThresMCPrompt_Cone_%d_Pt%d",icone,ipt);
497 snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
498 fhPtThresIsolatedPrompt[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
499 fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
500 outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
501
502 snprintf(name, buffersize,"hPtFracMCPrompt_Cone_%d_Pt%d",icone,ipt);
503 snprintf(title, buffersize,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
504 fhPtFracIsolatedPrompt[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
505 fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
506 outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
507
508 snprintf(name, buffersize,"hPtThresMCFragmentation_Cone_%d_Pt%d",icone,ipt);
509 snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
510 fhPtThresIsolatedFragmentation[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
511 fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
512 outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
513
514 snprintf(name, buffersize,"hPtFracMCFragmentation_Cone_%d_Pt%d",icone,ipt);
515 snprintf(title, buffersize,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
516 fhPtFracIsolatedFragmentation[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
517 fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
518 outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
519
520 snprintf(name, buffersize,"hPtThresMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
521 snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
522 fhPtThresIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
523 fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
524 outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
525
526 snprintf(name, buffersize,"hPtFracMCPi0Decay_Cone_%d_Pt%d",icone,ipt);
527 snprintf(title, buffersize,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
528 fhPtFracIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
529 fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
530 outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
531
532 snprintf(name, buffersize,"hPtThresMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
533 snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
534 fhPtThresIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
535 fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
536 outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
537
538 snprintf(name, buffersize,"hPtFracMCOtherDecay_Cone_%d_Pt%d",icone,ipt);
539 snprintf(title, buffersize,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
540 fhPtFracIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
541 fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
542 outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
543
544 snprintf(name, buffersize,"hPtThresMCConversion_Cone_%d_Pt%d",icone,ipt);
545 snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
546 fhPtThresIsolatedConversion[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
547 fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
548 outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
549
550 snprintf(name, buffersize,"hPtFracMCConversion_Cone_%d_Pt%d",icone,ipt);
551 snprintf(title, buffersize,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
552 fhPtFracIsolatedConversion[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
553 fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
554 outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
555
556 snprintf(name, buffersize,"hPtThresMCUnknown_Cone_%d_Pt%d",icone,ipt);
557 snprintf(title, buffersize,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
558 fhPtThresIsolatedUnknown[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
559 fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
560 outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
561
562 snprintf(name, buffersize,"hPtFracMCUnknown_Cone_%d_Pt%d",icone,ipt);
563 snprintf(title, buffersize,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
564 fhPtFracIsolatedUnknown[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
565 fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
566 outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
567
568 }//Histos with MC
569
570 }//icone loop
571 }//ipt loop
572 }
573
574 //Keep neutral meson selection histograms if requiered
575 //Setting done in AliNeutralMesonSelection
576 if(fMakeInvMass && GetNeutralMesonSelection()){
577 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
578 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
579 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
580 delete nmsHistos;
581 }
582
583 return outputContainer ;
584
585}
586
587//__________________________________________________________________
588void AliAnaParticleIsolation::MakeAnalysisFillAOD()
589{
590 //Do analysis and fill aods
591 //Search for the isolated photon in fCalorimeter with pt > GetMinPt()
592
593 if(!GetInputAODBranch()){
594 printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
595 abort();
596 }
597
598 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
599 printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
600 abort();
601 }
602
603 Int_t n = 0, nfrac = 0;
604 Bool_t isolated = kFALSE ;
605 Bool_t decay = kFALSE;
606 Float_t coneptsum = 0 ;
607 TObjArray * pl = 0x0; ;
608
609 //Select the calorimeter for candidate isolation with neutral particles
610 if(fCalorimeter == "PHOS")
611 pl = GetPHOSClusters();
612 else if (fCalorimeter == "EMCAL")
613 pl = GetEMCALClusters();
614
615 //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
616 Double_t ptLeading = 0. ;
617 Int_t idLeading = -1 ;
618 TLorentzVector mom ;
619 Int_t naod = GetInputAODBranch()->GetEntriesFast();
620 if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
621
622 for(Int_t iaod = 0; iaod < naod; iaod++){
623 AliAODPWG4ParticleCorrelation * aodinput = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
624
625 //If too small or too large pt, skip
626 if(aodinput->Pt() < GetMinPt() || aodinput->Pt() > GetMaxPt() ) continue ;
627
628 //check if it is low pt trigger particle, then adjust the isolation method
629 if(aodinput->Pt() < GetIsolationCut()->GetPtThreshold() || aodinput->Pt() < GetIsolationCut()->GetSumPtThreshold())
630 continue ; //trigger should not come from underlying event
631
632 //vertex cut in case of mixing
633 if (GetMixedEvent()) {
634 Int_t evt=-1;
635 Int_t id =-1;
636 if (aodinput->GetCaloLabel(0) >= 0 ){
637 id=aodinput->GetCaloLabel(0);
638 if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
639 }
640 else if(aodinput->GetTrackLabel(0) >= 0 ){
641 id=aodinput->GetTrackLabel(0);
642 if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
643 }
644 else continue;
645
646 if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
647 return ;
648 }
649
650 //find the leading particles with highest momentum
651 if ((aodinput->Pt())>ptLeading) {
652 ptLeading = aodinput->Pt() ;
226b95ba 653 idLeading = iaod ;
1a31a9ab 654 }
226b95ba 655 aodinput->SetLeadingParticle(kFALSE);
1a31a9ab 656 }//finish searching for leading trigger particle
657
658 // Check isolation of leading particle
659 if(idLeading < 0) return;
226b95ba 660
1a31a9ab 661 AliAODPWG4ParticleCorrelation * aodinput = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(idLeading));
226b95ba 662 aodinput->SetLeadingParticle(kTRUE);
1a31a9ab 663
664 //Check invariant mass, if pi0, tag decay.
665 if(fMakeInvMass) decay = CheckInvMass(idLeading,aodinput);
666
667 //After cuts, study isolation
668 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
669 GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,GetReader(), kTRUE, aodinput, GetAODObjArrayName(), n,nfrac,coneptsum, isolated);
670 aodinput->SetIsolated(isolated);
1a31a9ab 671 aodinput->SetTagged(decay);
672
673 if(GetDebug() > 1) {
674 if(isolated)printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
675 printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
676 }
677
678}
679
680//__________________________________________________________________
681void AliAnaParticleIsolation::MakeAnalysisFillHistograms()
682{
683 //Do analysis and fill histograms
684 Int_t n = 0, nfrac = 0;
685 Bool_t isolated = kFALSE ;
686 Float_t coneptsum = 0 ;
687 //Loop on stored AOD
688 Int_t naod = GetInputAODBranch()->GetEntriesFast();
689 if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
690
691 //Get vertex for photon momentum calculation
692 Double_t vertex[]={0,0,0} ; //vertex ;
693 //Double_t vertex2[]={0,0,0} ; //vertex ;
694 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
695 GetReader()->GetVertex(vertex);
696 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
697 }
698
699 for(Int_t iaod = 0; iaod < naod ; iaod++){
700 AliAODPWG4ParticleCorrelation* aod = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
701
702 if(!aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
703
704 Bool_t isolation = aod->IsIsolated();
705 Bool_t dec = aod->IsTagged();
706 Float_t pt = aod->Pt();
707 Float_t phi = aod->Phi();
708 Float_t eta = aod->Eta();
709 Float_t conesize = GetIsolationCut()->GetConeSize();
710
711 //Recover reference arrays with clusters and tracks
712 TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
713 TObjArray * reftracks = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
714 //If too small or too large pt, skip
715 if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
716
717 // --- In case of redoing isolation from delta AOD ----
718 if(fMakeSeveralIC) {
719 //Analysis of multiple IC at same time
720 MakeSeveralICAnalysis(aod);
721 }
722 else if(fReMakeIC){
723 //In case a more strict IC is needed in the produced AOD
724 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
725 GetIsolationCut()->MakeIsolationCut(reftracks, refclusters, GetReader(), kFALSE, aod, "", n,nfrac,coneptsum, isolated);
726 fhConeSumPt->Fill(pt,coneptsum);
727 if(GetDebug() > 0) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Energy Sum in Isolation Cone %2.2f\n", coneptsum);
728 }
729 // --- -------------------------------------------- ----
730
731 //Fill pt distribution of particles in cone
732 //Tracks
733 coneptsum = 0;
734 Double_t sumptFR = 0. ;
735 TObjArray * trackList = GetCTSTracks() ;
736 for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++){
737 AliVTrack* track = (AliVTrack *) trackList->At(itrack);
738 //fill the histograms at forward range
739 if(!track){
740 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
741 continue;
742 }
743 Double_t dPhi = phi - track->Phi() + TMath::PiOver2();
744 Double_t dEta = eta - track->Eta();
745 Double_t arg = dPhi*dPhi + dEta*dEta;
746 if(TMath::Sqrt(arg) < conesize){
747 fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
748 sumptFR+=track->Pt();
749 }
750 dPhi = phi - track->Phi() - TMath::PiOver2();
751 arg = dPhi*dPhi + dEta*dEta;
752 if(TMath::Sqrt(arg) < conesize){
753 fhPtInFRCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
754 sumptFR+=track->Pt();
755 }
756 }
757 fhFRConeSumPt->Fill(pt,sumptFR);
758 if(reftracks){
759 for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++){
760 AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
761 fhPtInCone->Fill(pt,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
762 coneptsum+=track->Pt();
763 }
764 }
765
766 //CaloClusters
767 if(refclusters){
768 TLorentzVector mom ;
769 for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++){
770 AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
771 calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
772
773 fhPtInCone->Fill(pt, mom.Pt());
774 coneptsum+=mom.Pt();
775 }
776 }
777
778 if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsum);
779
780 if(!fReMakeIC) fhConeSumPt->Fill(pt,coneptsum);
781
782 if(isolation){
783
784 if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED, fill histograms\n", iaod);
785
786 fhPtIso ->Fill(pt);
787 fhPhiIso ->Fill(pt,phi);
788 fhEtaIso ->Fill(pt,eta);
789 if (dec) fhPtInvMassDecayIso->Fill(pt);
790
791 if(IsDataMC()){
792 Int_t tag =aod->GetTag();
793
794 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
795 fhPtIsoPrompt ->Fill(pt);
796 fhPhiIsoPrompt ->Fill(pt,phi);
797 fhEtaIsoPrompt ->Fill(pt,eta);
798 }
799 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
800 {
801 fhPtIsoFragmentation ->Fill(pt);
802 fhPhiIsoFragmentation ->Fill(pt,phi);
803 fhEtaIsoFragmentation ->Fill(pt,eta);
804 }
805 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
806 {
807 fhPtIsoPi0Decay ->Fill(pt);
808 fhPhiIsoPi0Decay ->Fill(pt,phi);
809 fhEtaIsoPi0Decay ->Fill(pt,eta);
810 }
811 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
812 {
813 fhPtIsoOtherDecay ->Fill(pt);
814 fhPhiIsoOtherDecay ->Fill(pt,phi);
815 fhEtaIsoOtherDecay ->Fill(pt,eta);
816 }
817 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
818 {
819 fhPtIsoConversion ->Fill(pt);
820 fhPhiIsoConversion ->Fill(pt,phi);
821 fhEtaIsoConversion ->Fill(pt,eta);
822 }
823 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
824 {
825 fhPtIsoMCPhoton ->Fill(pt);
826 }
827
828 else
829 {
830 fhPtIsoUnknown ->Fill(pt);
831 fhPhiIsoUnknown ->Fill(pt,phi);
832 fhEtaIsoUnknown ->Fill(pt,eta);
833 }
834 }//Histograms with MC
835
836 }//Isolated histograms
837
838 if(!isolation)
839 {
840 fhPtNoIso ->Fill(pt);
841 if (dec) fhPtInvMassDecayNoIso->Fill(pt);
842
843 if(IsDataMC()){
844 Int_t tag =aod->GetTag();
845 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
846 {
847 fhPtNoIsoPi0Decay->Fill(pt);
848 }
849 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
850 {
851 fhPtNoIsoPrompt->Fill(pt);
852 }
853 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
854 {
855 fhPtNoIsoMCPhoton->Fill(pt);
856 }
857
858 }
859 }
860
861 }// aod loop
862
863}
864
865//____________________________________________________________________________
866void AliAnaParticleIsolation::InitParameters()
867{
868
869 //Initialize the parameters of the analysis.
870 SetInputAODName("PWG4Particle");
871 SetAODObjArrayName("IsolationCone");
872 AddToHistogramsName("AnaIsolation_");
873
874 fCalorimeter = "PHOS" ;
875 fReMakeIC = kFALSE ;
876 fMakeSeveralIC = kFALSE ;
877 fMakeInvMass = kFALSE ;
878
879 //----------- Several IC-----------------
880 fNCones = 5 ;
881 fNPtThresFrac = 5 ;
882 fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
883 fPtThresholds[0] = 1.; fPtThresholds[1] = 2.; fPtThresholds[2] = 3.; fPtThresholds[3] = 4.; fPtThresholds[4] = 5.;
884 fPtFractions[0] = 0.05; fPtFractions[1] = 0.075; fPtFractions[2] = 0.1; fPtFractions[3] = 1.25; fPtFractions[4] = 1.5;
885
886 //------------- Histograms settings -------
887 fHistoNPtSumBins = 100 ;
888 fHistoPtSumMax = 50 ;
889 fHistoPtSumMin = 0. ;
890
891 fHistoNPtInConeBins = 100 ;
892 fHistoPtInConeMax = 50 ;
893 fHistoPtInConeMin = 0. ;
894
895}
896
897//__________________________________________________________________
898void AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph)
899{
900 //Isolation Cut Analysis for both methods and different pt cuts and cones
901 Float_t ptC = ph->Pt();
902 Int_t tag = ph->GetTag();
903
904 //Keep original setting used when filling AODs, reset at end of analysis
905 Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
906 Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
907 Float_t rorg = GetIsolationCut()->GetConeSize();
908
909 Float_t coneptsum = 0 ;
910 Int_t n[10][10];//[fNCones][fNPtThresFrac];
911 Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
912 Bool_t isolated = kFALSE;
913
914 //Loop on cone sizes
915 for(Int_t icone = 0; icone<fNCones; icone++){
916 GetIsolationCut()->SetConeSize(fConeSizes[icone]);
917 coneptsum = 0 ;
918
919 //Loop on ptthresholds
920 for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
921 n[icone][ipt]=0;
922 nfrac[icone][ipt]=0;
923 GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
924 GetIsolationCut()->MakeIsolationCut(ph->GetObjArray(GetAODObjArrayName()+"Tracks"),
925 ph->GetObjArray(GetAODObjArrayName()+"Clusters"),
926 GetReader(), kFALSE, ph, "",n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
927
928 //Normal ptThreshold cut
929 if(n[icone][ipt] == 0) {
930 fhPtThresIsolated[icone][ipt]->Fill(ptC);
931 if(IsDataMC()){
932 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;
933 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;
934 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
935 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;
936 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay)) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;
937 else fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC) ;
938 }
939 }
940
941 //Pt threshold on pt cand/ pt in cone fraction
942 if(nfrac[icone][ipt] == 0) {
943 fhPtFracIsolated[icone][ipt]->Fill(ptC);
944 if(IsDataMC()){
945 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;
946 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;
947 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;
948 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;
949 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay)) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;
950 else fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ;
951 }
952 }
953 }//pt thresh loop
954
955 //Sum in cone histograms
956 fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
957 if(IsDataMC()){
958 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
959 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
960 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
961 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
962 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay)) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
963 else fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
964 }
965
966 }//cone size loop
967
968 //Reset original parameters for AOD analysis
969 GetIsolationCut()->SetPtThreshold(ptthresorg);
970 GetIsolationCut()->SetPtFraction(ptfracorg);
971 GetIsolationCut()->SetConeSize(rorg);
972
973}
974
975//__________________________________________________________________
976void AliAnaParticleIsolation::Print(const Option_t * opt) const
977{
978
979 //Print some relevant parameters set for the analysis
980 if(! opt)
981 return;
982
983 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
984 AliAnaPartCorrBaseClass::Print(" ");
985
986 printf("ReMake Isolation = %d \n", fReMakeIC) ;
987 printf("Make Several Isolation = %d \n", fMakeSeveralIC) ;
988 printf("Calorimeter for isolation = %s \n", fCalorimeter.Data()) ;
989
990 if(fMakeSeveralIC){
991 printf("N Cone Sizes = %d\n", fNCones) ;
992 printf("Cone Sizes = \n") ;
993 for(Int_t i = 0; i < fNCones; i++)
994 printf(" %1.2f;", fConeSizes[i]) ;
995 printf(" \n") ;
996
997 printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
998 printf(" pT thresholds = \n") ;
999 for(Int_t i = 0; i < fNPtThresFrac; i++)
1000 printf(" %2.2f;", fPtThresholds[i]) ;
1001
1002 printf(" \n") ;
1003
1004 printf(" pT fractions = \n") ;
1005 for(Int_t i = 0; i < fNPtThresFrac; i++)
1006 printf(" %2.2f;", fPtFractions[i]) ;
1007
1008 }
1009
1010 printf("Histograms: %3.1f < pT sum < %3.1f, Nbin = %d\n", fHistoPtSumMin, fHistoPtSumMax, fHistoNPtSumBins);
1011 printf("Histograms: %3.1f < pT in cone < %3.1f, Nbin = %d\n", fHistoPtInConeMin, fHistoPtInConeMax, fHistoNPtInConeBins);
1012
1013 printf(" \n") ;
1014
1015}
1016