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