]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/AliAnaGammaDirect.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
appropriately modified builders for package and PAR library
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaDirect.cxx
CommitLineData
f9cea31c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id$ */
16
17/* History of cvs commits:
18 *
19 * $Log$
5e7c69e7 20 * Revision 1.9 2007/11/17 16:39:49 gustavo
21 * removed deleting of not owned data and deleting of histograms which are exported to the output file (MG)
22 *
3788def4 23 * Revision 1.8 2007/10/29 13:48:42 gustavo
24 * Corrected coding violations
25 *
4b707925 26 * Revision 1.6 2007/08/17 12:40:04 schutz
27 * New analysis classes by Gustavo Conesa
28 *
bdcfac30 29 * Revision 1.4.4.4 2007/07/26 10:32:09 schutz
30 * new analysis classes in the the new analysis framework
2a1d8a29 31 *
f9cea31c 32 *
33 */
34
35//_________________________________________________________________________
bdcfac30 36// Class for the prompt gamma analysis, isolation cut
f9cea31c 37//
38// Class created from old AliPHOSGammaJet
39// (see AliRoot versions previous Release 4-09)
40//
2c8efea8 41// -- Author: Gustavo Conesa (LNF-INFN)
f9cea31c 42//////////////////////////////////////////////////////////////////////////////
bdcfac30 43
44
45// --- ROOT system ---
f9cea31c 46#include <TParticle.h>
47#include <TH2.h>
bdcfac30 48#include <TList.h>
f9cea31c 49#include "Riostream.h"
4b707925 50
2c8efea8 51// --- Analysis system ---
4b707925 52#include "AliAnaGammaDirect.h"
f9cea31c 53#include "AliLog.h"
2c8efea8 54#include "AliCaloTrackReader.h"
55#include "AliIsolationCut.h"
56#include "AliNeutralMesonSelection.h"
4b707925 57
f9cea31c 58ClassImp(AliAnaGammaDirect)
bdcfac30 59
f9cea31c 60//____________________________________________________________________________
bdcfac30 61 AliAnaGammaDirect::AliAnaGammaDirect() :
2c8efea8 62 AliAnaBaseClass(), fDetector(""), fMakeIC(0), fReMakeIC(0),
63 fMakeSeveralIC(0), fMakeInvMass(0),
64 fhPtGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
65 //Several IC
66 fNCones(0),fNPtThresFrac(0), fConeSizes(), fPtThresholds(), fPtFractions(),
67 fhPtThresIsolated(), fhPtFracIsolated(), fhPtSumIsolated(),
68 //MC
69 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
70 fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(), fhPtSumIsolatedPrompt(),
71 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
72 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
73 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
74 fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(), fhPtSumIsolatedPi0Decay(),
75 fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
76 fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(), fhPtSumIsolatedOtherDecay(),
77 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
78 fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(), fhPtSumIsolatedConversion(),
79 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0),
80 fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(), fhPtSumIsolatedUnknown()
f9cea31c 81{
bdcfac30 82 //default ctor
83
463ee300 84 //Initialize parameters
85 InitParameters();
86
5a67f6fd 87 for(Int_t i = 0; i < 5 ; i++){
88 fConeSizes[i] = 0 ;
89 fhPtSumIsolated[i] = 0 ;
90
91 fhPtSumIsolatedPrompt[i] = 0 ;
92 fhPtSumIsolatedFragmentation[i] = 0 ;
93 fhPtSumIsolatedPi0Decay[i] = 0 ;
94 fhPtSumIsolatedOtherDecay[i] = 0 ;
95 fhPtSumIsolatedConversion[i] = 0 ;
96 fhPtSumIsolatedUnknown[i] = 0 ;
97
98 for(Int_t j = 0; j < 5 ; j++){
99 fhPtThresIsolated[i][j] = 0 ;
100 fhPtFracIsolated[i][j] = 0 ;
101
102 fhPtThresIsolatedPrompt[i][j] = 0 ;
103 fhPtThresIsolatedFragmentation[i][j] = 0 ;
104 fhPtThresIsolatedPi0Decay[i][j] = 0 ;
105 fhPtThresIsolatedOtherDecay[i][j] = 0 ;
106 fhPtThresIsolatedConversion[i][j] = 0 ;
107 fhPtThresIsolatedUnknown[i][j] = 0 ;
108
109 fhPtFracIsolatedPrompt[i][j] = 0 ;
110 fhPtFracIsolatedFragmentation[i][j] = 0 ;
111 fhPtFracIsolatedPi0Decay[i][j] = 0 ;
112 fhPtFracIsolatedOtherDecay[i][j] = 0 ;
113 fhPtFracIsolatedConversion[i][j] = 0 ;
114 fhPtFracIsolatedUnknown[i][j] = 0 ;
115
116 }
117 }
118
119 for(Int_t i = 0; i < 5 ; i++){
120 fPtFractions[i]= 0 ;
121 fPtThresholds[i]= 0 ;
122 }
123
124
f9cea31c 125}
126
f9cea31c 127//____________________________________________________________________________
128AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
2c8efea8 129 AliAnaBaseClass(g), fDetector(g.fDetector),
130 fMakeIC(g.fMakeIC), fReMakeIC(g.fReMakeIC),
131 fMakeSeveralIC(g.fMakeSeveralIC), fMakeInvMass(g.fMakeInvMass),
132 fhPtGamma(g.fhPtGamma),fhPhiGamma(g.fhPhiGamma),
133 fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
134 //Several IC
135 fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(), fPtFractions(),
136 fhPtThresIsolated(), fhPtFracIsolated(), fhPtSumIsolated(),
137 //MC
138 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
139 fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(), fhPtSumIsolatedPrompt(),
140 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
141 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
142 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
143 fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(), fhPtSumIsolatedPi0Decay(),
144 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
145 fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(), fhPtSumIsolatedOtherDecay(),
146 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
147 fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(), fhPtSumIsolatedConversion(),
148 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
149 fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(), fhPtSumIsolatedUnknown()
f9cea31c 150{
151 // cpy ctor
bdcfac30 152
2c8efea8 153 //Several IC
bdcfac30 154 for(Int_t i = 0; i < fNCones ; i++){
155 fConeSizes[i] = g.fConeSizes[i];
156 fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
2c8efea8 157
158 fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
159 fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
160 fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
161 fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
162 fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
163 fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
164
165 for(Int_t j = 0; j < fNPtThresFrac ; j++){
bdcfac30 166 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
2c8efea8 167 fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];
168
169 fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
170 fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
171 fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
172 fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
173 fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
174 fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
175
176 fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
177 fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
178 fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
179 fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
180 fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
181 fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
182
183 }
bdcfac30 184 }
185
2c8efea8 186 for(Int_t i = 0; i < fNPtThresFrac ; i++){
187 fPtFractions[i]= g.fPtFractions[i];
bdcfac30 188 fPtThresholds[i]= g.fPtThresholds[i];
2c8efea8 189 }
4b707925 190
f9cea31c 191}
192
463ee300 193//_________________________________________________________________________
2c8efea8 194AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & g)
463ee300 195{
196 // assignment operator
bdcfac30 197
2c8efea8 198 if(&g == this) return *this;
199
200 fMakeIC = g.fMakeIC ;
201 fReMakeIC = g.fReMakeIC ;
202 fMakeSeveralIC = g.fMakeSeveralIC ;
203 fMakeInvMass = g.fMakeInvMass ;
204 fDetector = g.fDetector ;
205
206 fhPtGamma = g.fhPtGamma ;
207 fhPhiGamma = g.fhPhiGamma ;
208 fhEtaGamma = g.fhEtaGamma ;
209 fhConeSumPt = g.fhConeSumPt ;
210
211 fhPtPrompt = g.fhPtPrompt;
212 fhPhiPrompt = g.fhPhiPrompt;
213 fhEtaPrompt = g.fhEtaPrompt;
214 fhPtFragmentation = g.fhPtFragmentation;
215 fhPhiFragmentation = g.fhPhiFragmentation;
216 fhEtaFragmentation = g.fhEtaFragmentation;
217 fhPtPi0Decay = g.fhPtPi0Decay;
218 fhPhiPi0Decay = g.fhPhiPi0Decay;
219 fhEtaPi0Decay = g.fhEtaPi0Decay;
220 fhPtOtherDecay = g.fhPtOtherDecay;
221 fhPhiOtherDecay = g.fhPhiOtherDecay;
222 fhEtaOtherDecay = g.fhEtaOtherDecay;
223 fhPtConversion = g. fhPtConversion;
224 fhPhiConversion = g.fhPhiConversion;
225 fhEtaConversion = g.fhEtaConversion;
226 fhPtUnknown = g.fhPtUnknown;
227 fhPhiUnknown = g.fhPhiUnknown;
228 fhEtaUnknown = g.fhEtaUnknown;
229
230 //Several IC
231 fNCones = g.fNCones ;
232 fNPtThresFrac = g.fNPtThresFrac ;
bdcfac30 233
234 for(Int_t i = 0; i < fNCones ; i++){
2c8efea8 235 fConeSizes[i] = g.fConeSizes[i];
236 fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
237
238 fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
239 fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
240 fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
241 fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
242 fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
243 fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
244
245 for(Int_t j = 0; j < fNPtThresFrac ; j++){
246 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;
247 fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;
248
249 fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
250 fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
251 fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
252 fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
253 fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
254 fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
255
256 fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
257 fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
258 fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
259 fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
260 fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
261 fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
262
263 }
bdcfac30 264 }
265
2c8efea8 266 for(Int_t i = 0; i < fNPtThresFrac ; i++){
267 fPtThresholds[i]= g.fPtThresholds[i];
268 fPtFractions[i]= g.fPtFractions[i];
269 }
4b707925 270
463ee300 271 return *this;
bdcfac30 272
463ee300 273}
274
f9cea31c 275//____________________________________________________________________________
276AliAnaGammaDirect::~AliAnaGammaDirect()
277{
2c8efea8 278 //dtor
279 //do not delete histograms
280
281 delete [] fConeSizes ;
282 delete [] fPtThresholds ;
283 delete [] fPtFractions ;
284
285}
286//_________________________________________________________________________
287Bool_t AliAnaGammaDirect::CheckInvMass(const Int_t icalo,const TLorentzVector mom, Double_t *vertex, TClonesArray * pl){
288 //Search if there is a companion decay photon to the candidate
289 // and discard it in such case
290 TLorentzVector mom2 ;
291 for(Int_t jcalo = 0; jcalo < pl->GetEntries(); jcalo++){
292 if(icalo==jcalo) continue ;
293 AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (pl->At(jcalo));
294
295 //Cluster selection, not charged, with photon id and in fidutial cut, fill TLorentz
296 if(!SelectCluster(calo, vertex, mom2)) continue ;
297
298 //Select good pair (good phit, pt cuts, aperture and invariant mass)
299 if(GetNeutralMesonSelection()->SelectPair(mom, mom2)){
300 if(GetDebug()>1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom+mom2).Pt(), (mom+mom2).Phi(), (mom+mom2).Eta());
301 return kTRUE ;
302 }
303 }//loop
304
305 return kFALSE;
4b707925 306
2c8efea8 307}
308//_________________________________________________________________________
309Int_t AliAnaGammaDirect::CheckOrigin(const Int_t label){
310 //Play with the MC stack if available
311 //Check origin of the candidates, good for PYTHIA
312
313 AliStack * stack = GetMCStack() ;
314 if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
315
316 if(label >= 0 && label < stack->GetNtrack()){
317 //Mother
318 TParticle * mom = stack->Particle(label);
319 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
320 Int_t mStatus = mom->GetStatusCode() ;
321 Int_t iParent = mom->GetFirstMother() ;
322 if(label < 8 ) printf("Mother is parton %d\n",iParent);
323
324 //GrandParent
325 TParticle * parent = new TParticle ;
326 Int_t pPdg = -1;
327 Int_t pStatus =-1;
328 if(iParent > 0){
329 parent = stack->Particle(iParent);
330 pPdg = TMath::Abs(parent->GetPdgCode());
331 pStatus = parent->GetStatusCode();
332 }
333 else
334 printf("Parent with label %d\n",iParent);
335
336 //return tag
337 if(mPdg == 22){
338 if(mStatus == 1){
339 if(iParent < 8) {
340 if(pPdg == 22) return kPrompt;
341 else return kFragmentation;
342 }
343 else if(pStatus == 11){
344 if(pPdg == 111) return kPi0Decay ;
345 else if (pPdg == 321) return kEtaDecay ;
346 else return kOtherDecay ;
347 }
348 }//Status 1 : Pythia generated
349 else if(mStatus == 0){
350 if(pPdg ==22 || pPdg ==11) return kConversion ;
351 if(pPdg == 111) return kPi0Decay ;
352 else if (pPdg == 221) return kEtaDecay ;
353 else return kOtherDecay ;
354 }//status 0 : geant generated
355 }//Mother gamma
356 else if(mPdg == 111) return kPi0 ;
357 else if(mPdg == 221) return kEta ;
358 else if(mPdg ==11){
359 if(mStatus == 0) return kConversion ;
360 else return kElectron ;
361 }
362 else return kUnknown;
363 }//Good label value
364 else{
365 if(label < 0 ) printf("*** bad label or no stack ***: label %d \n", label);
366 if(label >= stack->GetNtrack()) printf("*** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
367 return kUnknown;
368 }//Bad label
369
370 return kUnknown;
371
2a1d8a29 372}
373
374//________________________________________________________________________
bdcfac30 375TList * AliAnaGammaDirect::GetCreateOutputObjects()
2a1d8a29 376{
463ee300 377 // Create histograms to be saved in output file and
bdcfac30 378 // store them in outputContainer
379 TList * outputContainer = new TList() ;
380 outputContainer->SetName("DirectGammaHistos") ;
2c8efea8 381
2a1d8a29 382 //Histograms of highest gamma identified in Event
2c8efea8 383 fhPtGamma = new TH1F("hPtGamma","Number of #gamma over calorimeter",240,0,120);
384 fhPtGamma->SetYTitle("N");
385 fhPtGamma->SetXTitle("p_{T #gamma}(GeV/c)");
386 outputContainer->Add(fhPtGamma) ;
2a1d8a29 387
388 fhPhiGamma = new TH2F
2c8efea8 389 ("hPhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
2a1d8a29 390 fhPhiGamma->SetYTitle("#phi");
391 fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
bdcfac30 392 outputContainer->Add(fhPhiGamma) ;
2a1d8a29 393
394 fhEtaGamma = new TH2F
2c8efea8 395 ("hEtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
2a1d8a29 396 fhEtaGamma->SetYTitle("#eta");
397 fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
bdcfac30 398 outputContainer->Add(fhEtaGamma) ;
399
2c8efea8 400 if(!fMakeSeveralIC){
401 fhConeSumPt = new TH2F
402 ("hConePtSum","#Sigma p_{T} in cone ",200,0,120,100,0,100);
403 fhConeSumPt->SetYTitle("#Sigma p_{T}");
404 fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
405 outputContainer->Add(fhConeSumPt) ;
406 }
4b707925 407
2c8efea8 408 if(IsDataMC()){
409
410 fhPtPrompt = new TH1F("hPtPrompt","Number of #gamma over calorimeter",240,0,120);
411 fhPtPrompt->SetYTitle("N");
412 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
413 outputContainer->Add(fhPtPrompt) ;
414
415 fhPhiPrompt = new TH2F
416 ("hPhiPrompt","#phi_{#gamma}",200,0,120,200,0,7);
417 fhPhiPrompt->SetYTitle("#phi");
418 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
419 outputContainer->Add(fhPhiPrompt) ;
420
421 fhEtaPrompt = new TH2F
422 ("hEtaPrompt","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
423 fhEtaPrompt->SetYTitle("#eta");
424 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
425 outputContainer->Add(fhEtaPrompt) ;
426
427 fhPtFragmentation = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",240,0,120);
428 fhPtFragmentation->SetYTitle("N");
429 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
430 outputContainer->Add(fhPtFragmentation) ;
431
432 fhPhiFragmentation = new TH2F
433 ("hPhiFragmentation","#phi_{#gamma}",200,0,120,200,0,7);
434 fhPhiFragmentation->SetYTitle("#phi");
435 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
436 outputContainer->Add(fhPhiFragmentation) ;
437
438 fhEtaFragmentation = new TH2F
439 ("hEtaFragmentation","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
440 fhEtaFragmentation->SetYTitle("#eta");
441 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
442 outputContainer->Add(fhEtaFragmentation) ;
443
444 fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",240,0,120);
445 fhPtPi0Decay->SetYTitle("N");
446 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
447 outputContainer->Add(fhPtPi0Decay) ;
448
449 fhPhiPi0Decay = new TH2F
450 ("hPhiPi0Decay","#phi_{#gamma}",200,0,120,200,0,7);
451 fhPhiPi0Decay->SetYTitle("#phi");
452 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
453 outputContainer->Add(fhPhiPi0Decay) ;
454
455 fhEtaPi0Decay = new TH2F
456 ("hEtaPi0Decay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
457 fhEtaPi0Decay->SetYTitle("#eta");
458 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
459 outputContainer->Add(fhEtaPi0Decay) ;
460
461 fhPtOtherDecay = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",240,0,120);
462 fhPtOtherDecay->SetYTitle("N");
463 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
464 outputContainer->Add(fhPtOtherDecay) ;
465
466 fhPhiOtherDecay = new TH2F
467 ("hPhiOtherDecay","#phi_{#gamma}",200,0,120,200,0,7);
468 fhPhiOtherDecay->SetYTitle("#phi");
469 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
470 outputContainer->Add(fhPhiOtherDecay) ;
471
472 fhEtaOtherDecay = new TH2F
473 ("hEtaOtherDecay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
474 fhEtaOtherDecay->SetYTitle("#eta");
475 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
476 outputContainer->Add(fhEtaOtherDecay) ;
477
478 fhPtConversion = new TH1F("hPtConversion","Number of #gamma over calorimeter",240,0,120);
479 fhPtConversion->SetYTitle("N");
480 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
481 outputContainer->Add(fhPtConversion) ;
482
483 fhPhiConversion = new TH2F
484 ("hPhiConversion","#phi_{#gamma}",200,0,120,200,0,7);
485 fhPhiConversion->SetYTitle("#phi");
486 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
487 outputContainer->Add(fhPhiConversion) ;
488
489 fhEtaConversion = new TH2F
490 ("hEtaConversion","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
491 fhEtaConversion->SetYTitle("#eta");
492 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
493 outputContainer->Add(fhEtaConversion) ;
494
495 fhPtUnknown = new TH1F("hPtUnknown","Number of #gamma over calorimeter",240,0,120);
496 fhPtUnknown->SetYTitle("N");
497 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
498 outputContainer->Add(fhPtUnknown) ;
499
500 fhPhiUnknown = new TH2F
501 ("hPhiUnknown","#phi_{#gamma}",200,0,120,200,0,7);
502 fhPhiUnknown->SetYTitle("#phi");
503 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
504 outputContainer->Add(fhPhiUnknown) ;
505
506 fhEtaUnknown = new TH2F
507 ("hEtaUnknown","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
508 fhEtaUnknown->SetYTitle("#eta");
509 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
510 outputContainer->Add(fhEtaUnknown) ;
511 }//Histos with MC
512
513 if(fMakeSeveralIC){
bdcfac30 514 char name[128];
515 char title[128];
516 for(Int_t icone = 0; icone<fNCones; icone++){
2c8efea8 517 sprintf(name,"hPtSumIsolated_Cone_%d",icone);
bdcfac30 518 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
519 fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
520 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
521 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
522 outputContainer->Add(fhPtSumIsolated[icone]) ;
2c8efea8 523
524 if(IsDataMC()){
525 sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);
526 sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
527 fhPtSumIsolatedPrompt[icone] = new TH2F(name, title,240,0,120,120,0,10);
528 fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
529 fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
530 outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
531
532 sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);
533 sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
534 fhPtSumIsolatedFragmentation[icone] = new TH2F(name, title,240,0,120,120,0,10);
535 fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
536 fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
537 outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
538
539 sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);
540 sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
541 fhPtSumIsolatedPi0Decay[icone] = new TH2F(name, title,240,0,120,120,0,10);
542 fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
543 fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
544 outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
545
546 sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);
547 sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
548 fhPtSumIsolatedOtherDecay[icone] = new TH2F(name, title,240,0,120,120,0,10);
549 fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
550 fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
551 outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
552
553 sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);
554 sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
555 fhPtSumIsolatedConversion[icone] = new TH2F(name, title,240,0,120,120,0,10);
556 fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
557 fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
558 outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
559
560 sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);
561 sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
562 fhPtSumIsolatedUnknown[icone] = new TH2F(name, title,240,0,120,120,0,10);
563 fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
564 fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
565 outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
566
567 }//Histos with MC
568
569 for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){
570 sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);
bdcfac30 571 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
572 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
573 fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
574 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
2c8efea8 575
576 sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);
577 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
578 fhPtFracIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
579 fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
580 outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
581
582 if(IsDataMC()){
583 sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);
584 sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
585 fhPtThresIsolatedPrompt[icone][ipt] = new TH1F(name, title,240,0,120);
586 fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
587 outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
588
589 sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);
590 sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
591 fhPtFracIsolatedPrompt[icone][ipt] = new TH1F(name, title,240,0,120);
592 fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
593 outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
594
595 sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
596 sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
597 fhPtThresIsolatedFragmentation[icone][ipt] = new TH1F(name, title,240,0,120);
598 fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
599 outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
600
601 sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
602 sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
603 fhPtFracIsolatedFragmentation[icone][ipt] = new TH1F(name, title,240,0,120);
604 fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
605 outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
606
607 sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
608 sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
609 fhPtThresIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,240,0,120);
610 fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
611 outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
612
613 sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
614 sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
615 fhPtFracIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,240,0,120);
616 fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
617 outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
618
619 sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
620 sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
621 fhPtThresIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,240,0,120);
622 fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
623 outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
624
625 sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
626 sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
627 fhPtFracIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,240,0,120);
628 fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
629 outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
630
631 sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);
632 sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
633 fhPtThresIsolatedConversion[icone][ipt] = new TH1F(name, title,240,0,120);
634 fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
635 outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
636
637 sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);
638 sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
639 fhPtFracIsolatedConversion[icone][ipt] = new TH1F(name, title,240,0,120);
640 fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
641 outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
642
643 sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);
644 sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
645 fhPtThresIsolatedUnknown[icone][ipt] = new TH1F(name, title,240,0,120);
646 fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
647 outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
648
649 sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);
650 sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
651 fhPtFracIsolatedUnknown[icone][ipt] = new TH1F(name, title,240,0,120);
652 fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
653 outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
654
655 }//Histos with MC
656
bdcfac30 657 }//icone loop
658 }//ipt loop
f9cea31c 659 }
660
2c8efea8 661 //Keep neutral meson selection histograms if requiered
662 //Setting done in AliNeutralMesonSelection
663 if(fMakeInvMass && GetNeutralMesonSelection()){
664 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
665 cout<<"NMSHistos "<< nmsHistos<<endl;
666 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
667 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
668 }
4b707925 669
bdcfac30 670 return outputContainer ;
463ee300 671
f9cea31c 672}
673
2c8efea8 674//__________________________________________________________________
675void AliAnaGammaDirect::MakeAnalysisFillAOD()
f9cea31c 676{
2c8efea8 677 //Do analysis and fill aods
678 //Search for the isolated photon in fDetector with pt > GetMinPt()
f9cea31c 679
2c8efea8 680 //Fill CaloClusters if working with ESDs
681 //if(GetReader()->GetDataType() == AliCaloTrackReader::kESD) ConnectAODCaloClusters();
4b707925 682
2c8efea8 683 Int_t n = 0, nfrac = 0;
684 Bool_t isolated = kFALSE ;
685 Float_t coneptsum = 0 ;
686 TClonesArray * pl = new TClonesArray;
687
688 //Get vertex for photon momentum calculation
689 Double_t vertex[]={0,0,0} ; //vertex ;
690 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
691
692 //Select the detector of the photon
693 if(fDetector == "PHOS")
694 pl = GetAODPHOS();
695 else if (fDetector == "EMCAL")
696 pl = GetAODEMCAL();
697 //cout<<"Number of entries "<<pl->GetEntries()<<endl;
f9cea31c 698
2c8efea8 699 //Fill AODCaloClusters and AODParticleCorrelation with PHOS aods
700 TLorentzVector mom ;
701 for(Int_t icalo = 0; icalo < pl->GetEntries(); icalo++){
702 AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (pl->At(icalo));
703
704 //Cluster selection, not charged, with photon id and in fidutial cut
705 if(!SelectCluster(calo,vertex,mom)) continue ;
706
707 //If too small pt, skip
708 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
709
710 //Play with the MC stack if available
711 Int_t tag =-1;
712 if(IsDataMC()){
713 //Check origin of the candidates
714 tag = CheckOrigin(calo->GetLabel(0));
715 if(GetDebug() > 0) printf("Origin of candidate %d\n",tag);
716 }//Work with stack also
717
718 //Check invariant mass
719 Bool_t decay = kFALSE ;
720 if(fMakeInvMass) decay = CheckInvMass(icalo,mom,vertex,pl);
721 if(decay) continue ;
722
723 //Create AOD for analysis
724 AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
725 ph.SetLabel(calo->GetLabel(0));
726 ph.SetPdg(AliCaloPID::kPhoton);
727 ph.SetDetector(fDetector);
728 ph.SetTag(tag);
729 if(fMakeIC){
730 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
731 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)GetAODCTS(), (TSeqCollection*)pl,
732 vertex, kTRUE, &ph,icalo,-1,
733 n,nfrac,coneptsum, isolated);
734 if(isolated){
735 //cout<<"Isolated : E "<<mom.E()<<" ; Phi"<<mom.Phi()<< " ; Eta "<<mom.Eta()<<endl;
736 AddAODParticleCorrelation(ph);
737 }
4b707925 738 }
2c8efea8 739 else //Fill all if not isolation requested
740 AddAODParticleCorrelation(ph);
4b707925 741
2c8efea8 742 }//loop
743
744 if(GetDebug() > 1) printf("End fill AODs ");
745
f9cea31c 746}
747
2c8efea8 748//__________________________________________________________________
749void AliAnaGammaDirect::MakeAnalysisFillHistograms()
750{
751 //Do analysis and fill histograms
752 Int_t n = 0, nfrac = 0;
753 Bool_t isolated = kFALSE ;
754 Float_t coneptsum = 0 ;
755 //cout<<"MakeAnalysisFillHistograms"<<endl;
756
757 //Get vertex for photon momentum calculation
758 Double_t v[]={0,0,0} ; //vertex ;
759 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(v);
760
761 //Loop on stored AOD photons
762 Int_t naod = GetAODBranch()->GetEntriesFast();
763 if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
764 for(Int_t iaod = 0; iaod < naod ; iaod++){
765 AliAODParticleCorrelation* ph = dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
766 Int_t pdg = ph->GetPdg();
767
768 //Only isolate photons in detector fDetector
769 if(pdg != AliCaloPID::kPhoton || ph->GetDetector() != fDetector) continue;
770
771 if(fMakeSeveralIC) {
772 //Analysis of multiple IC at same time
773 MakeSeveralICAnalysis(ph,v);
774 continue;
775 }
776 else if(fReMakeIC){
777 //In case a more strict IC is needed in the produced AOD
778 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
779 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
780 (TSeqCollection*)ph->GetRefIsolationConeClusters(),
781 v, kFALSE, ph,-1,-1,
782 n,nfrac,coneptsum, isolated);
783 }
784
785 //Fill histograms (normal case)
786 if(!fReMakeIC || (fReMakeIC && isolated)){
787
788 //printf("Isolated Gamma: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta()) ;
789
790 //Fill prompt gamma histograms
791 Float_t ptcluster = ph->Pt();
792 Float_t phicluster = ph->Phi();
793 Float_t etacluster = ph->Eta();
794
795 fhPtGamma ->Fill(ptcluster);
796 fhPhiGamma ->Fill(ptcluster,phicluster);
797 fhEtaGamma ->Fill(ptcluster,etacluster);
798 fhConeSumPt->Fill(ptcluster,coneptsum);
799
800 if(IsDataMC()){
801 Int_t tag =ph->GetTag();
802
803 if(tag == kPrompt){
804 fhPtPrompt ->Fill(ptcluster);
805 fhPhiPrompt ->Fill(ptcluster,phicluster);
806 fhEtaPrompt ->Fill(ptcluster,etacluster);
807 }
808 else if(tag==kFragmentation)
809 {
810 fhPtFragmentation ->Fill(ptcluster);
811 fhPhiFragmentation ->Fill(ptcluster,phicluster);
812 fhEtaFragmentation ->Fill(ptcluster,etacluster);
813 }
814 else if(tag==kPi0Decay)
815 {
816 fhPtPi0Decay ->Fill(ptcluster);
817 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
818 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
819 }
820 else if(tag==kEtaDecay || tag==kOtherDecay)
821 {
822 fhPtOtherDecay ->Fill(ptcluster);
823 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
824 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
825 }
826 else if(tag==kConversion)
827 {
828 fhPtConversion ->Fill(ptcluster);
829 fhPhiConversion ->Fill(ptcluster,phicluster);
830 fhEtaConversion ->Fill(ptcluster,etacluster);
831 }
832 else{
833
834 fhPtUnknown ->Fill(ptcluster);
835 fhPhiUnknown ->Fill(ptcluster,phicluster);
836 fhEtaUnknown ->Fill(ptcluster,etacluster);
837 }
838 }//Histograms with MC
839
840 }
841 }// aod loop
842
843}
844
845//____________________________________________________________________________
2a1d8a29 846void AliAnaGammaDirect::InitParameters()
f9cea31c 847{
2c8efea8 848
2a1d8a29 849 //Initialize the parameters of the analysis.
f9cea31c 850
2c8efea8 851 fDetector = "PHOS" ;
852 fMakeIC = kTRUE;
853 fReMakeIC = kFALSE ;
854 fMakeSeveralIC = kFALSE ;
855 fMakeInvMass = kFALSE ;
f9cea31c 856
2c8efea8 857 //----------- Several IC-----------------
4b707925 858 fNCones = 5 ;
2c8efea8 859 fNPtThresFrac = 6 ;
4b707925 860 fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
3bb2c538 861 fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
2c8efea8 862 fPtFractions[0]=0.05; fPtFractions[1]=0.075; fPtFractions[2]=0.1; fPtFractions[3]=1.25; fPtFractions[4]=1.5;fPtFractions[5]=2.;
bdcfac30 863
f9cea31c 864}
865
866//__________________________________________________________________
2c8efea8 867void AliAnaGammaDirect::MakeSeveralICAnalysis(AliAODParticleCorrelation* ph, Double_t v[3])
868{
869 //Isolation Cut Analysis for both methods and different pt cuts and cones
870 Float_t ptC = ph->Pt();
871 Float_t phiC = ph->Phi();
872 Float_t etaC = ph->Eta();
873
874 fhPtGamma->Fill(ptC);
875 fhPhiGamma->Fill(ptC,ph->Phi());
876 fhEtaGamma->Fill(ptC,ph->Eta());
877 Int_t tag =ph->GetTag();
878
879 if(IsDataMC()){
880 if(tag == kPrompt){
881 fhPtPrompt ->Fill(ptC);
882 fhPhiPrompt ->Fill(ptC,phiC);
883 fhEtaPrompt ->Fill(ptC,etaC);
f9cea31c 884 }
2c8efea8 885 else if(tag==kFragmentation)
886 {
887 fhPtFragmentation ->Fill(ptC);
888 fhPhiFragmentation ->Fill(ptC,phiC);
889 fhEtaFragmentation ->Fill(ptC,etaC);
890 }
891 else if(tag==kPi0Decay)
892 {
893 fhPtPi0Decay ->Fill(ptC);
894 fhPhiPi0Decay ->Fill(ptC,phiC);
895 fhEtaPi0Decay ->Fill(ptC,etaC);
896 }
897 else if(tag==kEtaDecay || tag==kOtherDecay)
898 {
899 fhPtOtherDecay ->Fill(ptC);
900 fhPhiOtherDecay ->Fill(ptC,phiC);
901 fhEtaOtherDecay ->Fill(ptC,etaC);
902 }
903 else if(tag==kConversion)
904 {
905 fhPtConversion ->Fill(ptC);
906 fhPhiConversion ->Fill(ptC,phiC);
907 fhEtaConversion ->Fill(ptC,etaC);
908 }
909 else{
f9cea31c 910
2c8efea8 911 fhPtUnknown ->Fill(ptC);
912 fhPhiUnknown ->Fill(ptC,phiC);
913 fhEtaUnknown ->Fill(ptC,etaC);
914 }
915 }//Histograms with MC
916 //Keep original setting used when filling AODs, reset at end of analysis
917 Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
918 Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
919 Float_t rorg = GetIsolationCut()->GetConeSize();
920
921 Float_t coneptsum = 0 ;
922 Int_t n[10][10];//[fNCones][fNPtThresFrac];
923 Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
924 Bool_t isolated = kFALSE;
925
926 for(Int_t icone = 0; icone<fNCones; icone++){
927 GetIsolationCut()->SetConeSize(fConeSizes[icone]);
928 coneptsum = 0 ;
929 for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
930 n[icone][ipt]=0;
931 nfrac[icone][ipt]=0;
932 GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
933 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
934 (TSeqCollection*)ph->GetRefIsolationConeClusters(),
935 v, kFALSE, ph,-1,-1,
936 n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
937 if(n[icone][ipt] == 0) {
938 fhPtThresIsolated[icone][ipt]->Fill(ptC);
939 if(IsDataMC()){
940 if(tag==kPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
941 else if(tag==kConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
942 else if(tag==kFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
943 else if(tag==kPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
944 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
945 else fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
946 }
947 }
948 if(nfrac[icone][ipt] == 0) {
949 fhPtFracIsolated[icone][ipt]->Fill(ptC);
950 if(IsDataMC()){
951 if(tag==kPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
952 else if(tag==kConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
953 else if(tag==kFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
954 else if(tag==kPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
955 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
956 else fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
957 }
f9cea31c 958 }
2c8efea8 959 }//pt thresh loop
960 fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
961 if(IsDataMC()){
962 if(tag==kPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
963 else if(tag==kConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
964 else if(tag==kFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
965 else if(tag==kPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
966 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
967 else fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
f9cea31c 968 }
f9cea31c 969
2c8efea8 970 }//cone size loop
971
972 //Reset original parameters for AOD analysis
973 GetIsolationCut()->SetPtThreshold(ptthresorg);
974 GetIsolationCut()->SetPtFraction(ptfracorg);
975 GetIsolationCut()->SetConeSize(rorg);
f9cea31c 976
bdcfac30 977}
978
4b707925 979//__________________________________________________________________
bdcfac30 980void AliAnaGammaDirect::Print(const Option_t * opt) const
981{
982
f9cea31c 983 //Print some relevant parameters set for the analysis
984 if(! opt)
985 return;
bdcfac30 986
2c8efea8 987 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
bdcfac30 988
2c8efea8 989 printf("Min Gamma pT = %2.2f\n", GetMinPt()) ;
990 printf("Max Gamma pT = %3.2f\n", GetMaxPt()) ;
991
992// if(IsCaloPIDOn())printf("Check PID \n") ;
993// if(IsCaloPIDRecalculationOn())printf("Recalculate PID \n") ;
994// if(IsFidutialCutOn())printf("Check Fidutial cut \n") ;
995 printf("Make Isolation = %d \n", fMakeIC) ;
996 printf("ReMake Isolation = %d \n", fReMakeIC) ;
997 printf("Make Several Isolation = %d \n", fMakeSeveralIC) ;
998
999 if(fMakeSeveralIC){
1000 printf("N Cone Sizes = %d\n", fNCones) ;
1001 printf("Cone Sizes = \n") ;
bdcfac30 1002 for(Int_t i = 0; i < fNCones; i++)
2c8efea8 1003 printf(" %1.2f;", fConeSizes[i]) ;
bdcfac30 1004 printf(" \n") ;
2c8efea8 1005
1006 printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
1007 printf(" pT thresholds = \n") ;
1008 for(Int_t i = 0; i < fNPtThresFrac; i++)
1009 printf(" %2.2f;", fPtThresholds[i]) ;
1010
1011 printf(" \n") ;
1012
1013 printf(" pT fractions = \n") ;
1014 for(Int_t i = 0; i < fNPtThresFrac; i++)
1015 printf(" %2.2f;", fPtFractions[i]) ;
1016
1017 }
1018
bdcfac30 1019 printf(" \n") ;
1020
1021}
2c8efea8 1022
1023//____________________________________________________________________________
1024Bool_t AliAnaGammaDirect::SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom){
1025 //Select cluster depending on its pid and acceptance selections
1026
1027 //Skip matched clusters with tracks
1028 if(calo->GetNTracksMatched() > 0) return kFALSE ;
1029
1030 //Check PID
1031 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1032 Int_t pdg = AliCaloPID::kPhoton;
1033 if(IsCaloPIDOn()){
1034 //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1035 //or redo PID, recommended option for EMCal.
1036 if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1037 pdg = GetCaloPID()->GetPdg(fDetector,calo->PID(),mom.E());//PID with weights
1038 else
1039 pdg = GetCaloPID()->GetPdg(fDetector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
1040
1041 if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
1042
1043 //If it does not pass pid, skip
1044 if(pdg != AliCaloPID::kPhoton) return kFALSE ;
1045 }
1046
1047 //Check acceptance selection
1048 if(IsFidutialCutOn()){
1049 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
1050 if(! in ) return kFALSE ;
1051 }
1052
1053 if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1054
1055 return kTRUE;
1056
1057 }