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 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 **************************************************************************/
17 /* History of cvs commits:
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)
23 * Revision 1.8 2007/10/29 13:48:42 gustavo
24 * Corrected coding violations
26 * Revision 1.6 2007/08/17 12:40:04 schutz
27 * New analysis classes by Gustavo Conesa
29 * Revision 1.4.4.4 2007/07/26 10:32:09 schutz
30 * new analysis classes in the the new analysis framework
35 //_________________________________________________________________________
36 // Class for the prompt gamma analysis, isolation cut
38 // Class created from old AliPHOSGammaJet
39 // (see AliRoot versions previous Release 4-09)
41 // -- Author: Gustavo Conesa (LNF-INFN)
42 //////////////////////////////////////////////////////////////////////////////
45 // --- ROOT system ---
46 #include <TParticle.h>
49 #include "Riostream.h"
51 // --- Analysis system ---
52 #include "AliAnaGammaDirect.h"
54 #include "AliCaloTrackReader.h"
55 #include "AliIsolationCut.h"
56 #include "AliNeutralMesonSelection.h"
58 ClassImp(AliAnaGammaDirect)
60 //____________________________________________________________________________
61 AliAnaGammaDirect::AliAnaGammaDirect() :
62 AliAnaBaseClass(), fDetector(""), fMakeIC(0), fReMakeIC(0),
63 fMakeSeveralIC(0), fMakeInvMass(0),
64 fhPtGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
66 fNCones(0),fNPtThresFrac(0), fConeSizes(), fPtThresholds(), fPtFractions(),
67 fhPtThresIsolated(), fhPtFracIsolated(), fhPtSumIsolated(),
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()
84 //Initialize parameters
89 //____________________________________________________________________________
90 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
91 AliAnaBaseClass(g), fDetector(g.fDetector),
92 fMakeIC(g.fMakeIC), fReMakeIC(g.fReMakeIC),
93 fMakeSeveralIC(g.fMakeSeveralIC), fMakeInvMass(g.fMakeInvMass),
94 fhPtGamma(g.fhPtGamma),fhPhiGamma(g.fhPhiGamma),
95 fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
97 fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(), fPtFractions(),
98 fhPtThresIsolated(), fhPtFracIsolated(), fhPtSumIsolated(),
100 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
101 fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(), fhPtSumIsolatedPrompt(),
102 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
103 fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(), fhPtSumIsolatedFragmentation(),
104 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
105 fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(), fhPtSumIsolatedPi0Decay(),
106 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
107 fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(), fhPtSumIsolatedOtherDecay(),
108 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
109 fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(), fhPtSumIsolatedConversion(),
110 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
111 fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(), fhPtSumIsolatedUnknown()
116 for(Int_t i = 0; i < fNCones ; i++){
117 fConeSizes[i] = g.fConeSizes[i];
118 fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
120 fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
121 fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
122 fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
123 fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
124 fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
125 fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
127 for(Int_t j = 0; j < fNPtThresFrac ; j++){
128 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
129 fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];
131 fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
132 fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
133 fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
134 fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
135 fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
136 fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
138 fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
139 fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
140 fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
141 fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
142 fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
143 fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
148 for(Int_t i = 0; i < fNPtThresFrac ; i++){
149 fPtFractions[i]= g.fPtFractions[i];
150 fPtThresholds[i]= g.fPtThresholds[i];
155 //_________________________________________________________________________
156 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & g)
158 // assignment operator
160 if(&g == this) return *this;
162 fMakeIC = g.fMakeIC ;
163 fReMakeIC = g.fReMakeIC ;
164 fMakeSeveralIC = g.fMakeSeveralIC ;
165 fMakeInvMass = g.fMakeInvMass ;
166 fDetector = g.fDetector ;
168 fhPtGamma = g.fhPtGamma ;
169 fhPhiGamma = g.fhPhiGamma ;
170 fhEtaGamma = g.fhEtaGamma ;
171 fhConeSumPt = g.fhConeSumPt ;
173 fhPtPrompt = g.fhPtPrompt;
174 fhPhiPrompt = g.fhPhiPrompt;
175 fhEtaPrompt = g.fhEtaPrompt;
176 fhPtFragmentation = g.fhPtFragmentation;
177 fhPhiFragmentation = g.fhPhiFragmentation;
178 fhEtaFragmentation = g.fhEtaFragmentation;
179 fhPtPi0Decay = g.fhPtPi0Decay;
180 fhPhiPi0Decay = g.fhPhiPi0Decay;
181 fhEtaPi0Decay = g.fhEtaPi0Decay;
182 fhPtOtherDecay = g.fhPtOtherDecay;
183 fhPhiOtherDecay = g.fhPhiOtherDecay;
184 fhEtaOtherDecay = g.fhEtaOtherDecay;
185 fhPtConversion = g. fhPtConversion;
186 fhPhiConversion = g.fhPhiConversion;
187 fhEtaConversion = g.fhEtaConversion;
188 fhPtUnknown = g.fhPtUnknown;
189 fhPhiUnknown = g.fhPhiUnknown;
190 fhEtaUnknown = g.fhEtaUnknown;
193 fNCones = g.fNCones ;
194 fNPtThresFrac = g.fNPtThresFrac ;
196 for(Int_t i = 0; i < fNCones ; i++){
197 fConeSizes[i] = g.fConeSizes[i];
198 fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
200 fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
201 fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
202 fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
203 fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
204 fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
205 fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
207 for(Int_t j = 0; j < fNPtThresFrac ; j++){
208 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;
209 fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;
211 fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
212 fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
213 fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
214 fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
215 fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
216 fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
218 fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
219 fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
220 fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
221 fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
222 fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
223 fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
228 for(Int_t i = 0; i < fNPtThresFrac ; i++){
229 fPtThresholds[i]= g.fPtThresholds[i];
230 fPtFractions[i]= g.fPtFractions[i];
237 //____________________________________________________________________________
238 AliAnaGammaDirect::~AliAnaGammaDirect()
241 //do not delete histograms
243 delete [] fConeSizes ;
244 delete [] fPtThresholds ;
245 delete [] fPtFractions ;
248 //_________________________________________________________________________
249 Bool_t AliAnaGammaDirect::CheckInvMass(const Int_t icalo,const TLorentzVector mom, Double_t *vertex, TClonesArray * pl){
250 //Search if there is a companion decay photon to the candidate
251 // and discard it in such case
252 TLorentzVector mom2 ;
253 for(Int_t jcalo = 0; jcalo < pl->GetEntries(); jcalo++){
254 if(icalo==jcalo) continue ;
255 AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (pl->At(jcalo));
257 //Cluster selection, not charged, with photon id and in fidutial cut, fill TLorentz
258 if(!SelectCluster(calo, vertex, mom2)) continue ;
260 //Select good pair (good phit, pt cuts, aperture and invariant mass)
261 if(GetNeutralMesonSelection()->SelectPair(mom, mom2)){
262 if(GetDebug()>1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom+mom2).Pt(), (mom+mom2).Phi(), (mom+mom2).Eta());
270 //_________________________________________________________________________
271 Int_t AliAnaGammaDirect::CheckOrigin(const Int_t label){
272 //Play with the MC stack if available
273 //Check origin of the candidates, good for PYTHIA
275 AliStack * stack = GetMCStack() ;
276 if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
278 if(label >= 0 && label < stack->GetNtrack()){
280 TParticle * mom = stack->Particle(label);
281 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
282 Int_t mStatus = mom->GetStatusCode() ;
283 Int_t iParent = mom->GetFirstMother() ;
284 if(label < 8 ) printf("Mother is parton %d\n",iParent);
287 TParticle * parent = new TParticle ;
291 parent = stack->Particle(iParent);
292 pPdg = TMath::Abs(parent->GetPdgCode());
293 pStatus = parent->GetStatusCode();
296 printf("Parent with label %d\n",iParent);
302 if(pPdg == 22) return kPrompt;
303 else return kFragmentation;
305 else if(pStatus == 11){
306 if(pPdg == 111) return kPi0Decay ;
307 else if (pPdg == 321) return kEtaDecay ;
308 else return kOtherDecay ;
310 }//Status 1 : Pythia generated
311 else if(mStatus == 0){
312 if(pPdg ==22 || pPdg ==11) return kConversion ;
313 if(pPdg == 111) return kPi0Decay ;
314 else if (pPdg == 221) return kEtaDecay ;
315 else return kOtherDecay ;
316 }//status 0 : geant generated
318 else if(mPdg == 111) return kPi0 ;
319 else if(mPdg == 221) return kEta ;
321 if(mStatus == 0) return kConversion ;
322 else return kElectron ;
324 else return kUnknown;
327 if(label < 0 ) printf("*** bad label or no stack ***: label %d \n", label);
328 if(label >= stack->GetNtrack()) printf("*** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
336 //________________________________________________________________________
337 TList * AliAnaGammaDirect::GetCreateOutputObjects()
339 // Create histograms to be saved in output file and
340 // store them in outputContainer
341 TList * outputContainer = new TList() ;
342 outputContainer->SetName("DirectGammaHistos") ;
344 //Histograms of highest gamma identified in Event
345 fhPtGamma = new TH1F("hPtGamma","Number of #gamma over calorimeter",240,0,120);
346 fhPtGamma->SetYTitle("N");
347 fhPtGamma->SetXTitle("p_{T #gamma}(GeV/c)");
348 outputContainer->Add(fhPtGamma) ;
350 fhPhiGamma = new TH2F
351 ("hPhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
352 fhPhiGamma->SetYTitle("#phi");
353 fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
354 outputContainer->Add(fhPhiGamma) ;
356 fhEtaGamma = new TH2F
357 ("hEtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
358 fhEtaGamma->SetYTitle("#eta");
359 fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
360 outputContainer->Add(fhEtaGamma) ;
363 fhConeSumPt = new TH2F
364 ("hConePtSum","#Sigma p_{T} in cone ",200,0,120,100,0,100);
365 fhConeSumPt->SetYTitle("#Sigma p_{T}");
366 fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
367 outputContainer->Add(fhConeSumPt) ;
372 fhPtPrompt = new TH1F("hPtPrompt","Number of #gamma over calorimeter",240,0,120);
373 fhPtPrompt->SetYTitle("N");
374 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
375 outputContainer->Add(fhPtPrompt) ;
377 fhPhiPrompt = new TH2F
378 ("hPhiPrompt","#phi_{#gamma}",200,0,120,200,0,7);
379 fhPhiPrompt->SetYTitle("#phi");
380 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
381 outputContainer->Add(fhPhiPrompt) ;
383 fhEtaPrompt = new TH2F
384 ("hEtaPrompt","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
385 fhEtaPrompt->SetYTitle("#eta");
386 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
387 outputContainer->Add(fhEtaPrompt) ;
389 fhPtFragmentation = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",240,0,120);
390 fhPtFragmentation->SetYTitle("N");
391 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
392 outputContainer->Add(fhPtFragmentation) ;
394 fhPhiFragmentation = new TH2F
395 ("hPhiFragmentation","#phi_{#gamma}",200,0,120,200,0,7);
396 fhPhiFragmentation->SetYTitle("#phi");
397 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
398 outputContainer->Add(fhPhiFragmentation) ;
400 fhEtaFragmentation = new TH2F
401 ("hEtaFragmentation","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
402 fhEtaFragmentation->SetYTitle("#eta");
403 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
404 outputContainer->Add(fhEtaFragmentation) ;
406 fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",240,0,120);
407 fhPtPi0Decay->SetYTitle("N");
408 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
409 outputContainer->Add(fhPtPi0Decay) ;
411 fhPhiPi0Decay = new TH2F
412 ("hPhiPi0Decay","#phi_{#gamma}",200,0,120,200,0,7);
413 fhPhiPi0Decay->SetYTitle("#phi");
414 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
415 outputContainer->Add(fhPhiPi0Decay) ;
417 fhEtaPi0Decay = new TH2F
418 ("hEtaPi0Decay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
419 fhEtaPi0Decay->SetYTitle("#eta");
420 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
421 outputContainer->Add(fhEtaPi0Decay) ;
423 fhPtOtherDecay = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",240,0,120);
424 fhPtOtherDecay->SetYTitle("N");
425 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
426 outputContainer->Add(fhPtOtherDecay) ;
428 fhPhiOtherDecay = new TH2F
429 ("hPhiOtherDecay","#phi_{#gamma}",200,0,120,200,0,7);
430 fhPhiOtherDecay->SetYTitle("#phi");
431 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
432 outputContainer->Add(fhPhiOtherDecay) ;
434 fhEtaOtherDecay = new TH2F
435 ("hEtaOtherDecay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
436 fhEtaOtherDecay->SetYTitle("#eta");
437 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
438 outputContainer->Add(fhEtaOtherDecay) ;
440 fhPtConversion = new TH1F("hPtConversion","Number of #gamma over calorimeter",240,0,120);
441 fhPtConversion->SetYTitle("N");
442 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
443 outputContainer->Add(fhPtConversion) ;
445 fhPhiConversion = new TH2F
446 ("hPhiConversion","#phi_{#gamma}",200,0,120,200,0,7);
447 fhPhiConversion->SetYTitle("#phi");
448 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
449 outputContainer->Add(fhPhiConversion) ;
451 fhEtaConversion = new TH2F
452 ("hEtaConversion","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
453 fhEtaConversion->SetYTitle("#eta");
454 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
455 outputContainer->Add(fhEtaConversion) ;
457 fhPtUnknown = new TH1F("hPtUnknown","Number of #gamma over calorimeter",240,0,120);
458 fhPtUnknown->SetYTitle("N");
459 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
460 outputContainer->Add(fhPtUnknown) ;
462 fhPhiUnknown = new TH2F
463 ("hPhiUnknown","#phi_{#gamma}",200,0,120,200,0,7);
464 fhPhiUnknown->SetYTitle("#phi");
465 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
466 outputContainer->Add(fhPhiUnknown) ;
468 fhEtaUnknown = new TH2F
469 ("hEtaUnknown","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
470 fhEtaUnknown->SetYTitle("#eta");
471 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
472 outputContainer->Add(fhEtaUnknown) ;
478 for(Int_t icone = 0; icone<fNCones; icone++){
479 sprintf(name,"hPtSumIsolated_Cone_%d",icone);
480 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
481 fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
482 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
483 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
484 outputContainer->Add(fhPtSumIsolated[icone]) ;
487 sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);
488 sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
489 fhPtSumIsolatedPrompt[icone] = new TH2F(name, title,240,0,120,120,0,10);
490 fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
491 fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
492 outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
494 sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);
495 sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
496 fhPtSumIsolatedFragmentation[icone] = new TH2F(name, title,240,0,120,120,0,10);
497 fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
498 fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
499 outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
501 sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);
502 sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
503 fhPtSumIsolatedPi0Decay[icone] = new TH2F(name, title,240,0,120,120,0,10);
504 fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
505 fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
506 outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
508 sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);
509 sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
510 fhPtSumIsolatedOtherDecay[icone] = new TH2F(name, title,240,0,120,120,0,10);
511 fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
512 fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
513 outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
515 sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);
516 sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
517 fhPtSumIsolatedConversion[icone] = new TH2F(name, title,240,0,120,120,0,10);
518 fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
519 fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
520 outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
522 sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);
523 sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
524 fhPtSumIsolatedUnknown[icone] = new TH2F(name, title,240,0,120,120,0,10);
525 fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
526 fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
527 outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
531 for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){
532 sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);
533 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
534 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
535 fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
536 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
538 sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);
539 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
540 fhPtFracIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
541 fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
542 outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
545 sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);
546 sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
547 fhPtThresIsolatedPrompt[icone][ipt] = new TH1F(name, title,240,0,120);
548 fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
549 outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
551 sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);
552 sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
553 fhPtFracIsolatedPrompt[icone][ipt] = new TH1F(name, title,240,0,120);
554 fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
555 outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
557 sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
558 sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
559 fhPtThresIsolatedFragmentation[icone][ipt] = new TH1F(name, title,240,0,120);
560 fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
561 outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
563 sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
564 sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
565 fhPtFracIsolatedFragmentation[icone][ipt] = new TH1F(name, title,240,0,120);
566 fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
567 outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
569 sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
570 sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
571 fhPtThresIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,240,0,120);
572 fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
573 outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
575 sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
576 sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
577 fhPtFracIsolatedPi0Decay[icone][ipt] = new TH1F(name, title,240,0,120);
578 fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
579 outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
581 sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
582 sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
583 fhPtThresIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,240,0,120);
584 fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
585 outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
587 sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
588 sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
589 fhPtFracIsolatedOtherDecay[icone][ipt] = new TH1F(name, title,240,0,120);
590 fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
591 outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
593 sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);
594 sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
595 fhPtThresIsolatedConversion[icone][ipt] = new TH1F(name, title,240,0,120);
596 fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
597 outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
599 sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);
600 sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
601 fhPtFracIsolatedConversion[icone][ipt] = new TH1F(name, title,240,0,120);
602 fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
603 outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
605 sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);
606 sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
607 fhPtThresIsolatedUnknown[icone][ipt] = new TH1F(name, title,240,0,120);
608 fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
609 outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
611 sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);
612 sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
613 fhPtFracIsolatedUnknown[icone][ipt] = new TH1F(name, title,240,0,120);
614 fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
615 outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
623 //Keep neutral meson selection histograms if requiered
624 //Setting done in AliNeutralMesonSelection
625 if(fMakeInvMass && GetNeutralMesonSelection()){
626 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
627 cout<<"NMSHistos "<< nmsHistos<<endl;
628 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
629 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
632 return outputContainer ;
636 //__________________________________________________________________
637 void AliAnaGammaDirect::MakeAnalysisFillAOD()
639 //Do analysis and fill aods
640 //Search for the isolated photon in fDetector with pt > GetMinPt()
642 //Fill CaloClusters if working with ESDs
643 //if(GetReader()->GetDataType() == AliCaloTrackReader::kESD) ConnectAODCaloClusters();
645 Int_t n = 0, nfrac = 0;
646 Bool_t isolated = kFALSE ;
647 Float_t coneptsum = 0 ;
648 TClonesArray * pl = new TClonesArray;
650 //Get vertex for photon momentum calculation
651 Double_t vertex[]={0,0,0} ; //vertex ;
652 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
654 //Select the detector of the photon
655 if(fDetector == "PHOS")
657 else if (fDetector == "EMCAL")
659 //cout<<"Number of entries "<<pl->GetEntries()<<endl;
661 //Fill AODCaloClusters and AODParticleCorrelation with PHOS aods
663 for(Int_t icalo = 0; icalo < pl->GetEntries(); icalo++){
664 AliAODCaloCluster * calo = dynamic_cast<AliAODCaloCluster*> (pl->At(icalo));
666 //Cluster selection, not charged, with photon id and in fidutial cut
667 if(!SelectCluster(calo,vertex,mom)) continue ;
669 //If too small pt, skip
670 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
672 //Play with the MC stack if available
675 //Check origin of the candidates
676 tag = CheckOrigin(calo->GetLabel(0));
677 if(GetDebug() > 0) printf("Origin of candidate %d\n",tag);
678 }//Work with stack also
680 //Check invariant mass
681 Bool_t decay = kFALSE ;
682 if(fMakeInvMass) decay = CheckInvMass(icalo,mom,vertex,pl);
685 //Create AOD for analysis
686 AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
687 ph.SetLabel(calo->GetLabel(0));
688 ph.SetPdg(AliCaloPID::kPhoton);
689 ph.SetDetector(fDetector);
692 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
693 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)GetAODCTS(), (TSeqCollection*)pl,
694 vertex, kTRUE, &ph,icalo,-1,
695 n,nfrac,coneptsum, isolated);
697 //cout<<"Isolated : E "<<mom.E()<<" ; Phi"<<mom.Phi()<< " ; Eta "<<mom.Eta()<<endl;
698 AddAODParticleCorrelation(ph);
701 else //Fill all if not isolation requested
702 AddAODParticleCorrelation(ph);
706 if(GetDebug() > 1) printf("End fill AODs ");
710 //__________________________________________________________________
711 void AliAnaGammaDirect::MakeAnalysisFillHistograms()
713 //Do analysis and fill histograms
714 Int_t n = 0, nfrac = 0;
715 Bool_t isolated = kFALSE ;
716 Float_t coneptsum = 0 ;
717 //cout<<"MakeAnalysisFillHistograms"<<endl;
719 //Get vertex for photon momentum calculation
720 Double_t v[]={0,0,0} ; //vertex ;
721 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(v);
723 //Loop on stored AOD photons
724 Int_t naod = GetAODBranch()->GetEntriesFast();
725 if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
726 for(Int_t iaod = 0; iaod < naod ; iaod++){
727 AliAODParticleCorrelation* ph = dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
728 Int_t pdg = ph->GetPdg();
730 //Only isolate photons in detector fDetector
731 if(pdg != AliCaloPID::kPhoton || ph->GetDetector() != fDetector) continue;
734 //Analysis of multiple IC at same time
735 MakeSeveralICAnalysis(ph,v);
739 //In case a more strict IC is needed in the produced AOD
740 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
741 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
742 (TSeqCollection*)ph->GetRefIsolationConeClusters(),
744 n,nfrac,coneptsum, isolated);
747 //Fill histograms (normal case)
748 if(!fReMakeIC || (fReMakeIC && isolated)){
750 //printf("Isolated Gamma: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta()) ;
752 //Fill prompt gamma histograms
753 Float_t ptcluster = ph->Pt();
754 Float_t phicluster = ph->Phi();
755 Float_t etacluster = ph->Eta();
757 fhPtGamma ->Fill(ptcluster);
758 fhPhiGamma ->Fill(ptcluster,phicluster);
759 fhEtaGamma ->Fill(ptcluster,etacluster);
760 fhConeSumPt->Fill(ptcluster,coneptsum);
763 Int_t tag =ph->GetTag();
766 fhPtPrompt ->Fill(ptcluster);
767 fhPhiPrompt ->Fill(ptcluster,phicluster);
768 fhEtaPrompt ->Fill(ptcluster,etacluster);
770 else if(tag==kFragmentation)
772 fhPtFragmentation ->Fill(ptcluster);
773 fhPhiFragmentation ->Fill(ptcluster,phicluster);
774 fhEtaFragmentation ->Fill(ptcluster,etacluster);
776 else if(tag==kPi0Decay)
778 fhPtPi0Decay ->Fill(ptcluster);
779 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
780 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
782 else if(tag==kEtaDecay || tag==kOtherDecay)
784 fhPtOtherDecay ->Fill(ptcluster);
785 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
786 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
788 else if(tag==kConversion)
790 fhPtConversion ->Fill(ptcluster);
791 fhPhiConversion ->Fill(ptcluster,phicluster);
792 fhEtaConversion ->Fill(ptcluster,etacluster);
796 fhPtUnknown ->Fill(ptcluster);
797 fhPhiUnknown ->Fill(ptcluster,phicluster);
798 fhEtaUnknown ->Fill(ptcluster,etacluster);
800 }//Histograms with MC
807 //____________________________________________________________________________
808 void AliAnaGammaDirect::InitParameters()
811 //Initialize the parameters of the analysis.
816 fMakeSeveralIC = kFALSE ;
817 fMakeInvMass = kFALSE ;
819 //----------- Several IC-----------------
822 fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
823 fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
824 fPtFractions[0]=0.05; fPtFractions[1]=0.075; fPtFractions[2]=0.1; fPtFractions[3]=1.25; fPtFractions[4]=1.5;fPtFractions[5]=2.;
828 //__________________________________________________________________
829 void AliAnaGammaDirect::MakeSeveralICAnalysis(AliAODParticleCorrelation* ph, Double_t v[3])
831 //Isolation Cut Analysis for both methods and different pt cuts and cones
832 Float_t ptC = ph->Pt();
833 Float_t phiC = ph->Phi();
834 Float_t etaC = ph->Eta();
836 fhPtGamma->Fill(ptC);
837 fhPhiGamma->Fill(ptC,ph->Phi());
838 fhEtaGamma->Fill(ptC,ph->Eta());
839 Int_t tag =ph->GetTag();
843 fhPtPrompt ->Fill(ptC);
844 fhPhiPrompt ->Fill(ptC,phiC);
845 fhEtaPrompt ->Fill(ptC,etaC);
847 else if(tag==kFragmentation)
849 fhPtFragmentation ->Fill(ptC);
850 fhPhiFragmentation ->Fill(ptC,phiC);
851 fhEtaFragmentation ->Fill(ptC,etaC);
853 else if(tag==kPi0Decay)
855 fhPtPi0Decay ->Fill(ptC);
856 fhPhiPi0Decay ->Fill(ptC,phiC);
857 fhEtaPi0Decay ->Fill(ptC,etaC);
859 else if(tag==kEtaDecay || tag==kOtherDecay)
861 fhPtOtherDecay ->Fill(ptC);
862 fhPhiOtherDecay ->Fill(ptC,phiC);
863 fhEtaOtherDecay ->Fill(ptC,etaC);
865 else if(tag==kConversion)
867 fhPtConversion ->Fill(ptC);
868 fhPhiConversion ->Fill(ptC,phiC);
869 fhEtaConversion ->Fill(ptC,etaC);
873 fhPtUnknown ->Fill(ptC);
874 fhPhiUnknown ->Fill(ptC,phiC);
875 fhEtaUnknown ->Fill(ptC,etaC);
877 }//Histograms with MC
878 //Keep original setting used when filling AODs, reset at end of analysis
879 Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
880 Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
881 Float_t rorg = GetIsolationCut()->GetConeSize();
883 Float_t coneptsum = 0 ;
884 Int_t n[10][10];//[fNCones][fNPtThresFrac];
885 Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
886 Bool_t isolated = kFALSE;
888 for(Int_t icone = 0; icone<fNCones; icone++){
889 GetIsolationCut()->SetConeSize(fConeSizes[icone]);
891 for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
894 GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
895 GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
896 (TSeqCollection*)ph->GetRefIsolationConeClusters(),
898 n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
899 if(n[icone][ipt] == 0) {
900 fhPtThresIsolated[icone][ipt]->Fill(ptC);
902 if(tag==kPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
903 else if(tag==kConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
904 else if(tag==kFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
905 else if(tag==kPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
906 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
907 else fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
910 if(nfrac[icone][ipt] == 0) {
911 fhPtFracIsolated[icone][ipt]->Fill(ptC);
913 if(tag==kPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
914 else if(tag==kConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
915 else if(tag==kFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
916 else if(tag==kPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
917 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
918 else fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
922 fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
924 if(tag==kPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
925 else if(tag==kConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
926 else if(tag==kFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
927 else if(tag==kPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
928 else if(tag==kOtherDecay || tag==kEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
929 else fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
934 //Reset original parameters for AOD analysis
935 GetIsolationCut()->SetPtThreshold(ptthresorg);
936 GetIsolationCut()->SetPtFraction(ptfracorg);
937 GetIsolationCut()->SetConeSize(rorg);
941 //__________________________________________________________________
942 void AliAnaGammaDirect::Print(const Option_t * opt) const
945 //Print some relevant parameters set for the analysis
949 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
951 printf("Min Gamma pT = %2.2f\n", GetMinPt()) ;
952 printf("Max Gamma pT = %3.2f\n", GetMaxPt()) ;
954 // if(IsCaloPIDOn())printf("Check PID \n") ;
955 // if(IsCaloPIDRecalculationOn())printf("Recalculate PID \n") ;
956 // if(IsFidutialCutOn())printf("Check Fidutial cut \n") ;
957 printf("Make Isolation = %d \n", fMakeIC) ;
958 printf("ReMake Isolation = %d \n", fReMakeIC) ;
959 printf("Make Several Isolation = %d \n", fMakeSeveralIC) ;
962 printf("N Cone Sizes = %d\n", fNCones) ;
963 printf("Cone Sizes = \n") ;
964 for(Int_t i = 0; i < fNCones; i++)
965 printf(" %1.2f;", fConeSizes[i]) ;
968 printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
969 printf(" pT thresholds = \n") ;
970 for(Int_t i = 0; i < fNPtThresFrac; i++)
971 printf(" %2.2f;", fPtThresholds[i]) ;
975 printf(" pT fractions = \n") ;
976 for(Int_t i = 0; i < fNPtThresFrac; i++)
977 printf(" %2.2f;", fPtFractions[i]) ;
985 //____________________________________________________________________________
986 Bool_t AliAnaGammaDirect::SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom){
987 //Select cluster depending on its pid and acceptance selections
989 //Skip matched clusters with tracks
990 if(calo->GetNTracksMatched() > 0) return kFALSE ;
993 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
994 Int_t pdg = AliCaloPID::kPhoton;
996 //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
997 //or redo PID, recommended option for EMCal.
998 if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
999 pdg = GetCaloPID()->GetPdg(fDetector,calo->PID(),mom.E());//PID with weights
1001 pdg = GetCaloPID()->GetPdg(fDetector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
1003 if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
1005 //If it does not pass pid, skip
1006 if(pdg != AliCaloPID::kPhoton) return kFALSE ;
1009 //Check acceptance selection
1010 if(IsFidutialCutOn()){
1011 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
1012 if(! in ) return kFALSE ;
1015 if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);