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