1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes iGetEntriesFast(s 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 **************************************************************************/
17 //_________________________________________________________________________
18 // Class for the prompt gamma analysis, isolation cut
20 // Class created from old AliPHOSGammaJet
21 // (see AliRoot versions previous Release 4-09)
23 // -- Author: Gustavo Conesa (LNF-INFN)
24 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
28 #include <TParticle.h>
31 #include "Riostream.h"
33 // --- Analysis system ---
34 #include "AliAnaGammaDirect.h"
36 #include "AliCaloTrackReader.h"
37 #include "AliIsolationCut.h"
38 #include "AliNeutralMesonSelection.h"
40 ClassImp(AliAnaGammaDirect)
42 //____________________________________________________________________________
43 AliAnaGammaDirect::AliAnaGammaDirect() :
44 AliAnaPartCorrBaseClass(), fDetector(""), fMakeIC(0), fReMakeIC(0),
45 fMakeSeveralIC(0), fMakeInvMass(0),
46 fhPtGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
48 fNCones(0),fNPtThresFrac(0), fConeSizes(), fPtThresholds(), fPtFractions(),
49 fhPtThresIsolated(), fhPtFracIsolated(), fhPtSumIsolated(),
51 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
52 fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(), fhPtSumIsolatedPrompt(),
53 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
54 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
55 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
56 fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(), fhPtSumIsolatedPi0Decay(),
57 fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
58 fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(), fhPtSumIsolatedOtherDecay(),
59 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
60 fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(), fhPtSumIsolatedConversion(),
61 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0),
62 fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(), fhPtSumIsolatedUnknown()
66 //Initialize parameters
69 for(Int_t i = 0; i < 5 ; i++){
71 fhPtSumIsolated[i] = 0 ;
73 fhPtSumIsolatedPrompt[i] = 0 ;
74 fhPtSumIsolatedFragmentation[i] = 0 ;
75 fhPtSumIsolatedPi0Decay[i] = 0 ;
76 fhPtSumIsolatedOtherDecay[i] = 0 ;
77 fhPtSumIsolatedConversion[i] = 0 ;
78 fhPtSumIsolatedUnknown[i] = 0 ;
80 for(Int_t j = 0; j < 5 ; j++){
81 fhPtThresIsolated[i][j] = 0 ;
82 fhPtFracIsolated[i][j] = 0 ;
84 fhPtThresIsolatedPrompt[i][j] = 0 ;
85 fhPtThresIsolatedFragmentation[i][j] = 0 ;
86 fhPtThresIsolatedPi0Decay[i][j] = 0 ;
87 fhPtThresIsolatedOtherDecay[i][j] = 0 ;
88 fhPtThresIsolatedConversion[i][j] = 0 ;
89 fhPtThresIsolatedUnknown[i][j] = 0 ;
91 fhPtFracIsolatedPrompt[i][j] = 0 ;
92 fhPtFracIsolatedFragmentation[i][j] = 0 ;
93 fhPtFracIsolatedPi0Decay[i][j] = 0 ;
94 fhPtFracIsolatedOtherDecay[i][j] = 0 ;
95 fhPtFracIsolatedConversion[i][j] = 0 ;
96 fhPtFracIsolatedUnknown[i][j] = 0 ;
101 for(Int_t i = 0; i < 5 ; i++){
103 fPtThresholds[i]= 0 ;
109 //____________________________________________________________________________
110 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
111 AliAnaPartCorrBaseClass(g), fDetector(g.fDetector),
112 fMakeIC(g.fMakeIC), fReMakeIC(g.fReMakeIC),
113 fMakeSeveralIC(g.fMakeSeveralIC), fMakeInvMass(g.fMakeInvMass),
114 fhPtGamma(g.fhPtGamma),fhPhiGamma(g.fhPhiGamma),
115 fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
117 fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(), fPtFractions(),
118 fhPtThresIsolated(), fhPtFracIsolated(), fhPtSumIsolated(),
120 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
121 fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(), fhPtSumIsolatedPrompt(),
122 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
123 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
124 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
125 fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(), fhPtSumIsolatedPi0Decay(),
126 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
127 fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(), fhPtSumIsolatedOtherDecay(),
128 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
129 fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(), fhPtSumIsolatedConversion(),
130 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
131 fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(), fhPtSumIsolatedUnknown()
136 for(Int_t i = 0; i < fNCones ; i++){
137 fConeSizes[i] = g.fConeSizes[i];
138 fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
140 fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
141 fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
142 fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
143 fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
144 fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
145 fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
147 for(Int_t j = 0; j < fNPtThresFrac ; j++){
148 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
149 fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];
151 fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
152 fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
153 fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
154 fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
155 fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
156 fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
158 fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
159 fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
160 fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
161 fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
162 fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
163 fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
168 for(Int_t i = 0; i < fNPtThresFrac ; i++){
169 fPtFractions[i]= g.fPtFractions[i];
170 fPtThresholds[i]= g.fPtThresholds[i];
175 //_________________________________________________________________________
176 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & g)
178 // assignment operator
180 if(&g == this) return *this;
182 fMakeIC = g.fMakeIC ;
183 fReMakeIC = g.fReMakeIC ;
184 fMakeSeveralIC = g.fMakeSeveralIC ;
185 fMakeInvMass = g.fMakeInvMass ;
186 fDetector = g.fDetector ;
188 fhPtGamma = g.fhPtGamma ;
189 fhPhiGamma = g.fhPhiGamma ;
190 fhEtaGamma = g.fhEtaGamma ;
191 fhConeSumPt = g.fhConeSumPt ;
193 fhPtPrompt = g.fhPtPrompt;
194 fhPhiPrompt = g.fhPhiPrompt;
195 fhEtaPrompt = g.fhEtaPrompt;
196 fhPtFragmentation = g.fhPtFragmentation;
197 fhPhiFragmentation = g.fhPhiFragmentation;
198 fhEtaFragmentation = g.fhEtaFragmentation;
199 fhPtPi0Decay = g.fhPtPi0Decay;
200 fhPhiPi0Decay = g.fhPhiPi0Decay;
201 fhEtaPi0Decay = g.fhEtaPi0Decay;
202 fhPtOtherDecay = g.fhPtOtherDecay;
203 fhPhiOtherDecay = g.fhPhiOtherDecay;
204 fhEtaOtherDecay = g.fhEtaOtherDecay;
205 fhPtConversion = g. fhPtConversion;
206 fhPhiConversion = g.fhPhiConversion;
207 fhEtaConversion = g.fhEtaConversion;
208 fhPtUnknown = g.fhPtUnknown;
209 fhPhiUnknown = g.fhPhiUnknown;
210 fhEtaUnknown = g.fhEtaUnknown;
213 fNCones = g.fNCones ;
214 fNPtThresFrac = g.fNPtThresFrac ;
216 for(Int_t i = 0; i < fNCones ; i++){
217 fConeSizes[i] = g.fConeSizes[i];
218 fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
220 fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
221 fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
222 fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
223 fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
224 fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
225 fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
227 for(Int_t j = 0; j < fNPtThresFrac ; j++){
228 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;
229 fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;
231 fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
232 fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
233 fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
234 fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
235 fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
236 fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
238 fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
239 fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
240 fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
241 fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
242 fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
243 fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
248 for(Int_t i = 0; i < fNPtThresFrac ; i++){
249 fPtThresholds[i]= g.fPtThresholds[i];
250 fPtFractions[i]= g.fPtFractions[i];
257 //____________________________________________________________________________
258 AliAnaGammaDirect::~AliAnaGammaDirect()
261 //do not delete histograms
263 delete [] fConeSizes ;
264 delete [] fPtThresholds ;
265 delete [] fPtFractions ;
268 //_________________________________________________________________________
269 Bool_t AliAnaGammaDirect::CheckInvMass(const Int_t icalo,const TLorentzVector mom, Double_t *vertex, TClonesArray * pl){
270 //Search if there is a companion decay photon to the candidate
271 // and discard it in such case
272 TLorentzVector mom2 ;
273 for(Int_t jcalo = 0; jcalo < pl->GetEntriesFast(); jcalo++){
274 if(icalo==jcalo) continue ;
275 AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (pl->At(jcalo));
277 //Cluster selection, not charged, with photon id and in fidutial cut, fill TLorentz
278 if(!SelectCluster(calo, vertex, mom2)) continue ;
280 //Select good pair (good phit, pt cuts, aperture and invariant mass)
281 if(GetNeutralMesonSelection()->SelectPair(mom, mom2)){
282 if(GetDebug()>1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom+mom2).Pt(), (mom+mom2).Phi(), (mom+mom2).Eta());
290 //_________________________________________________________________________
291 Int_t AliAnaGammaDirect::CheckOrigin(const Int_t label){
292 //Play with the MC stack if available
293 //Check origin of the candidates, good for PYTHIA
295 AliStack * stack = GetMCStack() ;
296 if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
298 if(label >= 0 && label < stack->GetNtrack()){
300 TParticle * mom = stack->Particle(label);
301 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
302 Int_t mStatus = mom->GetStatusCode() ;
303 Int_t iParent = mom->GetFirstMother() ;
304 if(label < 8 ) printf("Mother is parton %d\n",iParent);
307 TParticle * parent = new TParticle ;
311 parent = stack->Particle(iParent);
312 pPdg = TMath::Abs(parent->GetPdgCode());
313 pStatus = parent->GetStatusCode();
316 printf("Parent with label %d\n",iParent);
322 if(pPdg == 22) return kPrompt;
323 else return kFragmentation;
325 else if(pStatus == 11){
326 if(pPdg == 111) return kPi0Decay ;
327 else if (pPdg == 321) return kEtaDecay ;
328 else return kOtherDecay ;
330 }//Status 1 : Pythia generated
331 else if(mStatus == 0){
332 if(pPdg ==22 || pPdg ==11) return kConversion ;
333 if(pPdg == 111) return kPi0Decay ;
334 else if (pPdg == 221) return kEtaDecay ;
335 else return kOtherDecay ;
336 }//status 0 : geant generated
338 else if(mPdg == 111) return kPi0 ;
339 else if(mPdg == 221) return kEta ;
341 if(mStatus == 0) return kConversion ;
342 else return kElectron ;
344 else return kUnknown;
347 if(label < 0 ) printf("*** bad label or no stack ***: label %d \n", label);
348 if(label >= stack->GetNtrack()) printf("*** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
356 //________________________________________________________________________
357 TList * AliAnaGammaDirect::GetCreateOutputObjects()
359 // Create histograms to be saved in output file and
360 // store them in outputContainer
361 TList * outputContainer = new TList() ;
362 outputContainer->SetName("DirectGammaHistos") ;
364 //Histograms of highest gamma identified in Event
365 fhPtGamma = new TH1F("hPtGamma","Number of #gamma over calorimeter",240,0,120);
366 fhPtGamma->SetYTitle("N");
367 fhPtGamma->SetXTitle("p_{T #gamma}(GeV/c)");
368 outputContainer->Add(fhPtGamma) ;
370 fhPhiGamma = new TH2F
371 ("hPhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
372 fhPhiGamma->SetYTitle("#phi");
373 fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
374 outputContainer->Add(fhPhiGamma) ;
376 fhEtaGamma = new TH2F
377 ("hEtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
378 fhEtaGamma->SetYTitle("#eta");
379 fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
380 outputContainer->Add(fhEtaGamma) ;
383 fhConeSumPt = new TH2F
384 ("hConePtSum","#Sigma p_{T} in cone ",200,0,120,100,0,100);
385 fhConeSumPt->SetYTitle("#Sigma p_{T}");
386 fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
387 outputContainer->Add(fhConeSumPt) ;
392 fhPtPrompt = new TH1F("hPtPrompt","Number of #gamma over calorimeter",240,0,120);
393 fhPtPrompt->SetYTitle("N");
394 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
395 outputContainer->Add(fhPtPrompt) ;
397 fhPhiPrompt = new TH2F
398 ("hPhiPrompt","#phi_{#gamma}",200,0,120,200,0,7);
399 fhPhiPrompt->SetYTitle("#phi");
400 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
401 outputContainer->Add(fhPhiPrompt) ;
403 fhEtaPrompt = new TH2F
404 ("hEtaPrompt","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
405 fhEtaPrompt->SetYTitle("#eta");
406 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
407 outputContainer->Add(fhEtaPrompt) ;
409 fhPtFragmentation = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",240,0,120);
410 fhPtFragmentation->SetYTitle("N");
411 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
412 outputContainer->Add(fhPtFragmentation) ;
414 fhPhiFragmentation = new TH2F
415 ("hPhiFragmentation","#phi_{#gamma}",200,0,120,200,0,7);
416 fhPhiFragmentation->SetYTitle("#phi");
417 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
418 outputContainer->Add(fhPhiFragmentation) ;
420 fhEtaFragmentation = new TH2F
421 ("hEtaFragmentation","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
422 fhEtaFragmentation->SetYTitle("#eta");
423 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
424 outputContainer->Add(fhEtaFragmentation) ;
426 fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",240,0,120);
427 fhPtPi0Decay->SetYTitle("N");
428 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
429 outputContainer->Add(fhPtPi0Decay) ;
431 fhPhiPi0Decay = new TH2F
432 ("hPhiPi0Decay","#phi_{#gamma}",200,0,120,200,0,7);
433 fhPhiPi0Decay->SetYTitle("#phi");
434 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
435 outputContainer->Add(fhPhiPi0Decay) ;
437 fhEtaPi0Decay = new TH2F
438 ("hEtaPi0Decay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
439 fhEtaPi0Decay->SetYTitle("#eta");
440 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
441 outputContainer->Add(fhEtaPi0Decay) ;
443 fhPtOtherDecay = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",240,0,120);
444 fhPtOtherDecay->SetYTitle("N");
445 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
446 outputContainer->Add(fhPtOtherDecay) ;
448 fhPhiOtherDecay = new TH2F
449 ("hPhiOtherDecay","#phi_{#gamma}",200,0,120,200,0,7);
450 fhPhiOtherDecay->SetYTitle("#phi");
451 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
452 outputContainer->Add(fhPhiOtherDecay) ;
454 fhEtaOtherDecay = new TH2F
455 ("hEtaOtherDecay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
456 fhEtaOtherDecay->SetYTitle("#eta");
457 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
458 outputContainer->Add(fhEtaOtherDecay) ;
460 fhPtConversion = new TH1F("hPtConversion","Number of #gamma over calorimeter",240,0,120);
461 fhPtConversion->SetYTitle("N");
462 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
463 outputContainer->Add(fhPtConversion) ;
465 fhPhiConversion = new TH2F
466 ("hPhiConversion","#phi_{#gamma}",200,0,120,200,0,7);
467 fhPhiConversion->SetYTitle("#phi");
468 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
469 outputContainer->Add(fhPhiConversion) ;
471 fhEtaConversion = new TH2F
472 ("hEtaConversion","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
473 fhEtaConversion->SetYTitle("#eta");
474 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
475 outputContainer->Add(fhEtaConversion) ;
477 fhPtUnknown = new TH1F("hPtUnknown","Number of #gamma over calorimeter",240,0,120);
478 fhPtUnknown->SetYTitle("N");
479 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
480 outputContainer->Add(fhPtUnknown) ;
482 fhPhiUnknown = new TH2F
483 ("hPhiUnknown","#phi_{#gamma}",200,0,120,200,0,7);
484 fhPhiUnknown->SetYTitle("#phi");
485 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
486 outputContainer->Add(fhPhiUnknown) ;
488 fhEtaUnknown = new TH2F
489 ("hEtaUnknown","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
490 fhEtaUnknown->SetYTitle("#eta");
491 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
492 outputContainer->Add(fhEtaUnknown) ;
498 for(Int_t icone = 0; icone<fNCones; icone++){
499 sprintf(name,"hPtSumIsolated_Cone_%d",icone);
500 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
501 fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
502 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
503 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
504 outputContainer->Add(fhPtSumIsolated[icone]) ;
507 sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);
508 sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
509 fhPtSumIsolatedPrompt[icone] = new TH2F(name, title,240,0,120,120,0,10);
510 fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
511 fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
512 outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
514 sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);
515 sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
516 fhPtSumIsolatedFragmentation[icone] = new TH2F(name, title,240,0,120,120,0,10);
517 fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
518 fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
519 outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
521 sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);
522 sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
523 fhPtSumIsolatedPi0Decay[icone] = new TH2F(name, title,240,0,120,120,0,10);
524 fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
525 fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
526 outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
528 sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);
529 sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
530 fhPtSumIsolatedOtherDecay[icone] = new TH2F(name, title,240,0,120,120,0,10);
531 fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
532 fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
533 outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
535 sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);
536 sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
537 fhPtSumIsolatedConversion[icone] = new TH2F(name, title,240,0,120,120,0,10);
538 fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
539 fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
540 outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
542 sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);
543 sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
544 fhPtSumIsolatedUnknown[icone] = new TH2F(name, title,240,0,120,120,0,10);
545 fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
546 fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
547 outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
551 for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){
552 sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);
553 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
554 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
555 fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
556 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
558 sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);
559 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
560 fhPtFracIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
561 fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
562 outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
565 sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);
566 sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
567 fhPtThresIsolatedPrompt[icone][ipt] = new TH1F(name, title,240,0,120);
568 fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
569 outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
571 sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);
572 sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
573 fhPtFracIsolatedPrompt[icone][ipt] = new TH1F(name, title,240,0,120);
574 fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
575 outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
577 sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
578 sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
579 fhPtThresIsolatedFragmentation[icone][ipt] = new TH1F(name, title,240,0,120);
580 fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
581 outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
583 sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
584 sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
585 fhPtFracIsolatedFragmentation[icone][ipt] = new TH1F(name, title,240,0,120);
586 fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
587 outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
589 sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
590 sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
591 fhPtThresIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,240,0,120);
592 fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
593 outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
595 sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
596 sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
597 fhPtFracIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,240,0,120);
598 fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
599 outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
601 sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
602 sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
603 fhPtThresIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,240,0,120);
604 fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
605 outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
607 sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
608 sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
609 fhPtFracIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,240,0,120);
610 fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
611 outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
613 sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);
614 sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
615 fhPtThresIsolatedConversion[icone][ipt] = new TH1F(name, title,240,0,120);
616 fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
617 outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
619 sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);
620 sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
621 fhPtFracIsolatedConversion[icone][ipt] = new TH1F(name, title,240,0,120);
622 fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
623 outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
625 sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);
626 sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
627 fhPtThresIsolatedUnknown[icone][ipt] = new TH1F(name, title,240,0,120);
628 fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
629 outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
631 sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);
632 sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
633 fhPtFracIsolatedUnknown[icone][ipt] = new TH1F(name, title,240,0,120);
634 fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
635 outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
643 //Keep neutral meson selection histograms if requiered
644 //Setting done in AliNeutralMesonSelection
645 if(fMakeInvMass && GetNeutralMesonSelection()){
646 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
647 cout<<"NMSHistos "<< nmsHistos<<endl;
648 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
649 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
652 return outputContainer ;
656 //__________________________________________________________________
657 void AliAnaGammaDirect::MakeAnalysisFillAOD()
659 //Do analysis and fill aods
660 //Search for the isolated photon in fDetector with pt > GetMinPt()
662 //Fill CaloClusters if working with ESDs
663 //if(GetReader()->GetDataType() == AliCaloTrackReader::kESD) ConnectAODCaloClusters();
665 Int_t n = 0, nfrac = 0;
666 Bool_t isolated = kFALSE ;
667 Float_t coneptsum = 0 ;
668 TClonesArray * pl = new TClonesArray;
670 //Get vertex for photon momentum calculation
671 Double_t vertex[]={0,0,0} ; //vertex ;
672 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
674 //Select the detector of the photon
675 if(fDetector == "PHOS")
677 else if (fDetector == "EMCAL")
679 //cout<<"Number of entries "<<pl->GetEntriesFast()<<endl;
681 //Fill AODCaloClusters and AODParticleCorrelation with PHOS aods
683 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
684 AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (pl->At(icalo));
686 //Cluster selection, not charged, with photon id and in fidutial cut
687 if(!SelectCluster(calo,vertex,mom)) continue ;
689 //If too small pt, skip
690 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
692 //Play with the MC stack if available
695 //Check origin of the candidates
696 tag = CheckOrigin(calo->GetLabel(0));
697 if(GetDebug() > 0) printf("Origin of candidate %d\n",tag);
698 }//Work with stack also
700 //Check invariant mass
701 Bool_t decay = kFALSE ;
702 if(fMakeInvMass) decay = CheckInvMass(icalo,mom,vertex,pl);
705 //Create AOD for analysis
706 AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
707 ph.SetLabel(calo->GetLabel(0));
708 ph.SetPdg(AliCaloPID::kPhoton);
709 ph.SetDetector(fDetector);
712 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
713 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)GetAODCTS(), (TSeqCollection*)pl,
714 vertex, kTRUE, &ph,icalo,-1,
715 n,nfrac,coneptsum, isolated);
717 //cout<<"Isolated : E "<<mom.E()<<" ; Phi"<<mom.Phi()<< " ; Eta "<<mom.Eta()<<endl;
718 AddAODParticleCorrelation(ph);
721 else //Fill all if not isolation requested
722 AddAODParticleCorrelation(ph);
726 if(GetDebug() > 1) printf("End fill AODs ");
730 //__________________________________________________________________
731 void AliAnaGammaDirect::MakeAnalysisFillHistograms()
733 //Do analysis and fill histograms
734 Int_t n = 0, nfrac = 0;
735 Bool_t isolated = kFALSE ;
736 Float_t coneptsum = 0 ;
737 //cout<<"MakeAnalysisFillHistograms"<<endl;
739 //Get vertex for photon momentum calculation
740 Double_t v[]={0,0,0} ; //vertex ;
741 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(v);
743 //Loop on stored AOD photons
744 Int_t naod = GetAODBranch()->GetEntriesFast();
745 if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
746 for(Int_t iaod = 0; iaod < naod ; iaod++){
747 AliAODParticleCorrelation* ph = dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
748 Int_t pdg = ph->GetPdg();
750 //Only isolate photons in detector fDetector
751 if(pdg != AliCaloPID::kPhoton || ph->GetDetector() != fDetector) continue;
754 //Analysis of multiple IC at same time
755 MakeSeveralICAnalysis(ph,v);
759 //In case a more strict IC is needed in the produced AOD
760 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
761 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
762 (TSeqCollection*)ph->GetRefIsolationConeClusters(),
764 n,nfrac,coneptsum, isolated);
767 //Fill histograms (normal case)
768 if(!fReMakeIC || (fReMakeIC && isolated)){
770 //printf("Isolated Gamma: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta()) ;
772 //Fill prompt gamma histograms
773 Float_t ptcluster = ph->Pt();
774 Float_t phicluster = ph->Phi();
775 Float_t etacluster = ph->Eta();
777 fhPtGamma ->Fill(ptcluster);
778 fhPhiGamma ->Fill(ptcluster,phicluster);
779 fhEtaGamma ->Fill(ptcluster,etacluster);
780 fhConeSumPt->Fill(ptcluster,coneptsum);
783 Int_t tag =ph->GetTag();
786 fhPtPrompt ->Fill(ptcluster);
787 fhPhiPrompt ->Fill(ptcluster,phicluster);
788 fhEtaPrompt ->Fill(ptcluster,etacluster);
790 else if(tag==kFragmentation)
792 fhPtFragmentation ->Fill(ptcluster);
793 fhPhiFragmentation ->Fill(ptcluster,phicluster);
794 fhEtaFragmentation ->Fill(ptcluster,etacluster);
796 else if(tag==kPi0Decay)
798 fhPtPi0Decay ->Fill(ptcluster);
799 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
800 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
802 else if(tag==kEtaDecay || tag==kOtherDecay)
804 fhPtOtherDecay ->Fill(ptcluster);
805 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
806 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
808 else if(tag==kConversion)
810 fhPtConversion ->Fill(ptcluster);
811 fhPhiConversion ->Fill(ptcluster,phicluster);
812 fhEtaConversion ->Fill(ptcluster,etacluster);
816 fhPtUnknown ->Fill(ptcluster);
817 fhPhiUnknown ->Fill(ptcluster,phicluster);
818 fhEtaUnknown ->Fill(ptcluster,etacluster);
820 }//Histograms with MC
827 //____________________________________________________________________________
828 void AliAnaGammaDirect::InitParameters()
831 //Initialize the parameters of the analysis.
836 fMakeSeveralIC = kFALSE ;
837 fMakeInvMass = kFALSE ;
839 //----------- Several IC-----------------
842 fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
843 fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
844 fPtFractions[0]=0.05; fPtFractions[1]=0.075; fPtFractions[2]=0.1; fPtFractions[3]=1.25; fPtFractions[4]=1.5;fPtFractions[5]=2.;
848 //__________________________________________________________________
849 void AliAnaGammaDirect::MakeSeveralICAnalysis(AliAODParticleCorrelation* ph, Double_t v[3])
851 //Isolation Cut Analysis for both methods and different pt cuts and cones
852 Float_t ptC = ph->Pt();
853 Float_t phiC = ph->Phi();
854 Float_t etaC = ph->Eta();
856 fhPtGamma->Fill(ptC);
857 fhPhiGamma->Fill(ptC,ph->Phi());
858 fhEtaGamma->Fill(ptC,ph->Eta());
859 Int_t tag =ph->GetTag();
863 fhPtPrompt ->Fill(ptC);
864 fhPhiPrompt ->Fill(ptC,phiC);
865 fhEtaPrompt ->Fill(ptC,etaC);
867 else if(tag==kFragmentation)
869 fhPtFragmentation ->Fill(ptC);
870 fhPhiFragmentation ->Fill(ptC,phiC);
871 fhEtaFragmentation ->Fill(ptC,etaC);
873 else if(tag==kPi0Decay)
875 fhPtPi0Decay ->Fill(ptC);
876 fhPhiPi0Decay ->Fill(ptC,phiC);
877 fhEtaPi0Decay ->Fill(ptC,etaC);
879 else if(tag==kEtaDecay || tag==kOtherDecay)
881 fhPtOtherDecay ->Fill(ptC);
882 fhPhiOtherDecay ->Fill(ptC,phiC);
883 fhEtaOtherDecay ->Fill(ptC,etaC);
885 else if(tag==kConversion)
887 fhPtConversion ->Fill(ptC);
888 fhPhiConversion ->Fill(ptC,phiC);
889 fhEtaConversion ->Fill(ptC,etaC);
893 fhPtUnknown ->Fill(ptC);
894 fhPhiUnknown ->Fill(ptC,phiC);
895 fhEtaUnknown ->Fill(ptC,etaC);
897 }//Histograms with MC
898 //Keep original setting used when filling AODs, reset at end of analysis
899 Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
900 Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
901 Float_t rorg = GetIsolationCut()->GetConeSize();
903 Float_t coneptsum = 0 ;
904 Int_t n[10][10];//[fNCones][fNPtThresFrac];
905 Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
906 Bool_t isolated = kFALSE;
908 for(Int_t icone = 0; icone<fNCones; icone++){
909 GetIsolationCut()->SetConeSize(fConeSizes[icone]);
911 for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
914 GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
915 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
916 (TSeqCollection*)ph->GetRefIsolationConeClusters(),
918 n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
919 if(n[icone][ipt] == 0) {
920 fhPtThresIsolated[icone][ipt]->Fill(ptC);
922 if(tag==kPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
923 else if(tag==kConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
924 else if(tag==kFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
925 else if(tag==kPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
926 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
927 else fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
930 if(nfrac[icone][ipt] == 0) {
931 fhPtFracIsolated[icone][ipt]->Fill(ptC);
933 if(tag==kPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
934 else if(tag==kConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
935 else if(tag==kFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
936 else if(tag==kPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
937 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
938 else fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
942 fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
944 if(tag==kPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
945 else if(tag==kConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
946 else if(tag==kFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
947 else if(tag==kPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
948 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
949 else fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
954 //Reset original parameters for AOD analysis
955 GetIsolationCut()->SetPtThreshold(ptthresorg);
956 GetIsolationCut()->SetPtFraction(ptfracorg);
957 GetIsolationCut()->SetConeSize(rorg);
961 //__________________________________________________________________
962 void AliAnaGammaDirect::Print(const Option_t * opt) const
965 //Print some relevant parameters set for the analysis
969 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
971 printf("Min Gamma pT = %2.2f\n", GetMinPt()) ;
972 printf("Max Gamma pT = %3.2f\n", GetMaxPt()) ;
974 // if(IsCaloPIDOn())printf("Check PID \n") ;
975 // if(IsCaloPIDRecalculationOn())printf("Recalculate PID \n") ;
976 // if(IsFidutialCutOn())printf("Check Fidutial cut \n") ;
977 printf("Make Isolation = %d \n", fMakeIC) ;
978 printf("ReMake Isolation = %d \n", fReMakeIC) ;
979 printf("Make Several Isolation = %d \n", fMakeSeveralIC) ;
982 printf("N Cone Sizes = %d\n", fNCones) ;
983 printf("Cone Sizes = \n") ;
984 for(Int_t i = 0; i < fNCones; i++)
985 printf(" %1.2f;", fConeSizes[i]) ;
988 printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
989 printf(" pT thresholds = \n") ;
990 for(Int_t i = 0; i < fNPtThresFrac; i++)
991 printf(" %2.2f;", fPtThresholds[i]) ;
995 printf(" pT fractions = \n") ;
996 for(Int_t i = 0; i < fNPtThresFrac; i++)
997 printf(" %2.2f;", fPtFractions[i]) ;
1005 //____________________________________________________________________________
1006 Bool_t AliAnaGammaDirect::SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom){
1007 //Select cluster depending on its pid and acceptance selections
1009 //Skip matched clusters with tracks
1010 if(calo->GetNTracksMatched() > 0) return kFALSE ;
1013 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1014 Int_t pdg = AliCaloPID::kPhoton;
1016 //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1017 //or redo PID, recommended option for EMCal.
1018 if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1019 pdg = GetCaloPID()->GetPdg(fDetector,calo->PID(),mom.E());//PID with weights
1021 pdg = GetCaloPID()->GetPdg(fDetector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
1023 if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
1025 //If it does not pass pid, skip
1026 if(pdg != AliCaloPID::kPhoton) return kFALSE ;
1029 //Check acceptance selection
1030 if(IsFidutialCutOn()){
1031 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
1032 if(! in ) return kFALSE ;
1035 if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);