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