cf53dab96f4c01bfad8c691a2f5181e23485984f
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaDirect.cxx
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$
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  *
23  * Revision 1.8  2007/10/29 13:48:42  gustavo
24  * Corrected coding violations
25  *
26  * Revision 1.6  2007/08/17 12:40:04  schutz
27  * New analysis classes by Gustavo Conesa
28  *
29  * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
30  * new analysis classes in the the new analysis framework
31  *
32  *
33  */
34
35 //_________________________________________________________________________
36 // Class for the prompt gamma analysis, isolation cut
37 //
38 //  Class created from old AliPHOSGammaJet 
39 //  (see AliRoot versions previous Release 4-09)
40 //
41 // -- Author: Gustavo Conesa (LNF-INFN) 
42 //////////////////////////////////////////////////////////////////////////////
43   
44   
45 // --- ROOT system --- 
46 #include <TParticle.h>
47 #include <TH2.h>
48 #include <TList.h>
49 #include "Riostream.h"
50
51 // --- Analysis system --- 
52 #include "AliAnaGammaDirect.h" 
53 #include "AliLog.h"
54 #include "AliCaloTrackReader.h"
55 #include "AliIsolationCut.h"
56 #include "AliNeutralMesonSelection.h"
57
58 ClassImp(AliAnaGammaDirect)
59   
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),
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()
81 {
82   //default ctor
83   
84   //Initialize parameters
85   InitParameters();
86
87 }
88
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),  
96   //Several IC
97   fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(),  fPtFractions(), 
98   fhPtThresIsolated(),  fhPtFracIsolated(), fhPtSumIsolated(),
99   //MC
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()
112 {
113   // cpy ctor
114   
115   //Several IC
116   for(Int_t i = 0; i < fNCones ; i++){
117     fConeSizes[i] =  g.fConeSizes[i];
118     fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; 
119
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]; 
126
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];
130
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]; 
137  
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]; 
144
145     } 
146   }
147   
148   for(Int_t i = 0; i < fNPtThresFrac ; i++){
149     fPtFractions[i]=   g.fPtFractions[i];
150     fPtThresholds[i]=   g.fPtThresholds[i];
151   }
152
153 }
154
155 //_________________________________________________________________________
156 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & g)
157 {
158   // assignment operator
159   
160   if(&g == this) return *this;
161
162   fMakeIC = g.fMakeIC ;
163   fReMakeIC = g.fReMakeIC ;
164   fMakeSeveralIC = g.fMakeSeveralIC ;
165   fMakeInvMass = g.fMakeInvMass ;
166   fDetector = g.fDetector ;
167
168   fhPtGamma = g.fhPtGamma ; 
169   fhPhiGamma = g.fhPhiGamma ;
170   fhEtaGamma = g.fhEtaGamma ;
171   fhConeSumPt = g.fhConeSumPt ;
172
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; 
191
192   //Several IC
193   fNCones = g.fNCones ;
194   fNPtThresFrac = g.fNPtThresFrac ; 
195    
196   for(Int_t i = 0; i < fNCones ; i++){
197     fConeSizes[i] =  g.fConeSizes[i];
198     fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
199
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]; 
206
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] ;
210
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]; 
217  
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]; 
224
225     }
226   }
227   
228   for(Int_t i = 0; i < fNPtThresFrac ; i++){
229     fPtThresholds[i]=   g.fPtThresholds[i];
230     fPtFractions[i]=   g.fPtFractions[i];
231   }
232
233   return *this;
234   
235 }
236
237 //____________________________________________________________________________
238 AliAnaGammaDirect::~AliAnaGammaDirect() 
239 {
240   //dtor
241   //do not delete histograms
242
243  delete [] fConeSizes ; 
244  delete [] fPtThresholds ; 
245  delete [] fPtFractions ; 
246
247 }
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));
256   
257     //Cluster selection, not charged, with photon id and in fidutial cut, fill TLorentz
258     if(!SelectCluster(calo, vertex, mom2)) continue ;
259
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());
263       return kTRUE ;
264     }
265   }//loop
266
267   return kFALSE;
268
269 }
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
274   
275   AliStack * stack =  GetMCStack() ;
276   if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
277   
278   if(label >= 0 && label <  stack->GetNtrack()){
279     //Mother
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);
285     
286     //GrandParent
287     TParticle * parent = new TParticle ;
288     Int_t pPdg = -1;
289     Int_t pStatus =-1;
290     if(iParent > 0){
291       parent = stack->Particle(iParent);
292       pPdg = TMath::Abs(parent->GetPdgCode());
293       pStatus = parent->GetStatusCode();  
294     }
295     else
296       printf("Parent with label %d\n",iParent);
297     
298     //return tag
299     if(mPdg == 22){
300       if(mStatus == 1){
301         if(iParent < 8) {
302           if(pPdg == 22) return kPrompt;
303           else  return kFragmentation;
304         }
305         else if(pStatus == 11){
306           if(pPdg == 111) return kPi0Decay ;
307           else if (pPdg == 321)  return kEtaDecay ;
308           else  return kOtherDecay ;
309         }
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
317     }//Mother gamma
318     else if(mPdg == 111)  return kPi0 ;
319     else if(mPdg == 221)  return kEta ;
320     else if(mPdg ==11){
321       if(mStatus == 0) return kConversion ;
322       else return kElectron ;
323     }
324     else return kUnknown;
325   }//Good label value
326   else{
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());
329     return kUnknown;
330   }//Bad label
331   
332   return kUnknown;
333   
334 }
335
336 //________________________________________________________________________
337 TList *  AliAnaGammaDirect::GetCreateOutputObjects()
338 {  
339   // Create histograms to be saved in output file and 
340   // store them in outputContainer
341   TList * outputContainer = new TList() ; 
342   outputContainer->SetName("DirectGammaHistos") ; 
343
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) ; 
349   
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) ; 
355   
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) ;
361
362   if(!fMakeSeveralIC){
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) ;
368   }
369   
370   if(IsDataMC()){
371     
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) ; 
376     
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) ; 
382     
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) ;
388     
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) ; 
393     
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) ; 
399     
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) ;
405     
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) ; 
410     
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) ; 
416     
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) ;
422     
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) ; 
427     
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) ; 
433     
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) ;
439     
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) ; 
444     
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) ; 
450     
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) ;
456     
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) ; 
461     
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) ; 
467     
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) ;
473   }//Histos with MC
474   
475   if(fMakeSeveralIC){
476     char name[128];
477     char title[128];
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]) ; 
485     
486       if(IsDataMC()){
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]) ; 
493
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]) ; 
500
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]) ; 
507
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]) ; 
514
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]) ; 
521
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]) ; 
528
529       }//Histos with MC
530
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]) ; 
537
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]) ; 
543
544         if(IsDataMC()){
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]) ; 
550           
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]) ; 
556
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]) ; 
562           
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]) ; 
568
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]) ; 
574           
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]) ; 
580
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]) ; 
586           
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]) ;
592
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]) ; 
598           
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]) ;
604
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]) ; 
610           
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]) ;  
616  
617         }//Histos with MC
618
619       }//icone loop
620     }//ipt loop
621   }
622
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)) ;
630   }
631
632   return outputContainer ;
633
634 }
635
636 //__________________________________________________________________
637 void  AliAnaGammaDirect::MakeAnalysisFillAOD() 
638 {
639   //Do analysis and fill aods
640   //Search for the isolated photon in fDetector with pt > GetMinPt()
641
642   //Fill CaloClusters if working with ESDs
643   //if(GetReader()->GetDataType() == AliCaloTrackReader::kESD) ConnectAODCaloClusters(); 
644
645   Int_t n = 0, nfrac = 0;
646   Bool_t isolated = kFALSE ; 
647   Float_t coneptsum = 0 ;
648   TClonesArray * pl = new TClonesArray; 
649
650   //Get vertex for photon momentum calculation
651   Double_t vertex[]={0,0,0} ; //vertex ;
652   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
653
654   //Select the detector of the photon
655   if(fDetector == "PHOS")
656     pl = GetAODPHOS();
657   else if (fDetector == "EMCAL")
658     pl = GetAODEMCAL();
659   //cout<<"Number of entries "<<pl->GetEntries()<<endl;
660   
661   //Fill AODCaloClusters and AODParticleCorrelation with PHOS aods
662   TLorentzVector mom ;
663   for(Int_t icalo = 0; icalo < pl->GetEntries(); icalo++){
664     AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (pl->At(icalo));
665   
666     //Cluster selection, not charged, with photon id and in fidutial cut
667     if(!SelectCluster(calo,vertex,mom)) continue ;
668     
669     //If too small pt, skip
670     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
671
672     //Play with the MC stack if available
673     Int_t tag =-1;
674     if(IsDataMC()){
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   
679
680     //Check invariant mass
681     Bool_t decay = kFALSE ;
682     if(fMakeInvMass) decay = CheckInvMass(icalo,mom,vertex,pl);
683     if(decay) continue ;
684
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);
690     ph.SetTag(tag);  
691     if(fMakeIC){
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);
696       if(isolated){
697         //cout<<"Isolated : E "<<mom.E()<<" ; Phi"<<mom.Phi()<< " ; Eta "<<mom.Eta()<<endl;
698         AddAODParticleCorrelation(ph);
699       }
700     }
701     else //Fill all if not isolation requested
702       AddAODParticleCorrelation(ph);
703
704   }//loop
705   
706   if(GetDebug() > 1) printf("End fill AODs ");  
707   
708 }
709
710 //__________________________________________________________________
711 void  AliAnaGammaDirect::MakeAnalysisFillHistograms() 
712 {
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;
718
719   //Get vertex for photon momentum calculation
720   Double_t v[]={0,0,0} ; //vertex ;
721   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(v);
722
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();
729
730    //Only isolate photons in detector fDetector
731    if(pdg != AliCaloPID::kPhoton || ph->GetDetector() != fDetector) continue;
732    
733    if(fMakeSeveralIC) {
734      //Analysis of multiple IC at same time
735      MakeSeveralICAnalysis(ph,v);
736      continue;
737    }
738    else if(fReMakeIC){
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(), 
743                                          v, kFALSE, ph,-1,-1,
744                                          n,nfrac,coneptsum, isolated);
745    }
746
747    //Fill histograms (normal case)
748    if(!fReMakeIC || (fReMakeIC && isolated)){
749      
750      //printf("Isolated Gamma: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta()) ;
751      
752      //Fill prompt gamma histograms 
753      Float_t ptcluster = ph->Pt();
754      Float_t phicluster = ph->Phi();
755      Float_t etacluster = ph->Eta();
756      
757      fhPtGamma   ->Fill(ptcluster);
758      fhPhiGamma ->Fill(ptcluster,phicluster);
759      fhEtaGamma ->Fill(ptcluster,etacluster);
760      fhConeSumPt->Fill(ptcluster,coneptsum);
761
762      if(IsDataMC()){
763        Int_t tag =ph->GetTag();
764            
765        if(tag == kPrompt){
766          fhPtPrompt   ->Fill(ptcluster);
767          fhPhiPrompt ->Fill(ptcluster,phicluster);
768          fhEtaPrompt ->Fill(ptcluster,etacluster);
769        }
770        else if(tag==kFragmentation)
771          {
772            fhPtFragmentation   ->Fill(ptcluster);
773            fhPhiFragmentation ->Fill(ptcluster,phicluster);
774            fhEtaFragmentation ->Fill(ptcluster,etacluster);
775          }
776        else if(tag==kPi0Decay)
777          {
778            fhPtPi0Decay   ->Fill(ptcluster);
779            fhPhiPi0Decay ->Fill(ptcluster,phicluster);
780            fhEtaPi0Decay ->Fill(ptcluster,etacluster);
781          }
782        else if(tag==kEtaDecay || tag==kOtherDecay)
783          {
784            fhPtOtherDecay   ->Fill(ptcluster);
785            fhPhiOtherDecay ->Fill(ptcluster,phicluster);
786            fhEtaOtherDecay ->Fill(ptcluster,etacluster);
787          }
788        else if(tag==kConversion)
789          {
790            fhPtConversion   ->Fill(ptcluster);
791            fhPhiConversion ->Fill(ptcluster,phicluster);
792            fhEtaConversion ->Fill(ptcluster,etacluster);
793          }
794        else{
795          
796          fhPtUnknown   ->Fill(ptcluster);
797          fhPhiUnknown ->Fill(ptcluster,phicluster);
798          fhEtaUnknown ->Fill(ptcluster,etacluster);
799        }
800      }//Histograms with MC
801      
802    }
803   }// aod loop
804   
805 }
806
807 //____________________________________________________________________________
808 void AliAnaGammaDirect::InitParameters()
809 {
810   
811   //Initialize the parameters of the analysis.
812
813   fDetector = "PHOS" ;
814   fMakeIC = kTRUE;
815   fReMakeIC = kFALSE ;
816   fMakeSeveralIC = kFALSE ;
817   fMakeInvMass = kFALSE ;
818
819  //----------- Several IC-----------------
820   fNCones           = 5 ; 
821   fNPtThresFrac         = 6 ; 
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.;
825
826 }
827
828 //__________________________________________________________________
829 void  AliAnaGammaDirect::MakeSeveralICAnalysis(AliAODParticleCorrelation* ph, Double_t v[3]) 
830 {
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();
835
836   fhPtGamma->Fill(ptC);
837   fhPhiGamma->Fill(ptC,ph->Phi());
838   fhEtaGamma->Fill(ptC,ph->Eta());
839   Int_t tag =ph->GetTag();
840
841   if(IsDataMC()){
842     if(tag == kPrompt){
843       fhPtPrompt   ->Fill(ptC);
844       fhPhiPrompt ->Fill(ptC,phiC);
845       fhEtaPrompt ->Fill(ptC,etaC);
846     }
847     else if(tag==kFragmentation)
848       {
849         fhPtFragmentation   ->Fill(ptC);
850         fhPhiFragmentation ->Fill(ptC,phiC);
851         fhEtaFragmentation ->Fill(ptC,etaC);
852       }
853     else if(tag==kPi0Decay)
854       {
855         fhPtPi0Decay   ->Fill(ptC);
856         fhPhiPi0Decay ->Fill(ptC,phiC);
857         fhEtaPi0Decay ->Fill(ptC,etaC);
858       }
859     else if(tag==kEtaDecay || tag==kOtherDecay)
860       {
861         fhPtOtherDecay   ->Fill(ptC);
862         fhPhiOtherDecay ->Fill(ptC,phiC);
863         fhEtaOtherDecay ->Fill(ptC,etaC);
864       }
865     else if(tag==kConversion)
866       {
867         fhPtConversion   ->Fill(ptC);
868         fhPhiConversion ->Fill(ptC,phiC);
869         fhEtaConversion ->Fill(ptC,etaC);
870       }
871     else{
872       
873       fhPtUnknown   ->Fill(ptC);
874       fhPhiUnknown ->Fill(ptC,phiC);
875       fhEtaUnknown ->Fill(ptC,etaC);
876     }
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();
882
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;
887
888   for(Int_t icone = 0; icone<fNCones; icone++){
889     GetIsolationCut()->SetConeSize(fConeSizes[icone]);
890     coneptsum = 0 ;
891     for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
892       n[icone][ipt]=0;
893       nfrac[icone][ipt]=0;
894       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
895       GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(), 
896                                           (TSeqCollection*)ph->GetRefIsolationConeClusters(), 
897                                           v, kFALSE, ph,-1,-1,
898                                           n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
899       if(n[icone][ipt] == 0) {
900         fhPtThresIsolated[icone][ipt]->Fill(ptC);
901         if(IsDataMC()){
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) ;
908         }
909       }
910       if(nfrac[icone][ipt] == 0) {
911         fhPtFracIsolated[icone][ipt]->Fill(ptC);
912         if(IsDataMC()){
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) ;
919         }
920       }
921     }//pt thresh loop
922     fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
923     if(IsDataMC()){
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) ;
930     }
931
932   }//cone size loop
933
934   //Reset original parameters for AOD analysis
935   GetIsolationCut()->SetPtThreshold(ptthresorg);
936   GetIsolationCut()->SetPtFraction(ptfracorg);
937   GetIsolationCut()->SetConeSize(rorg);
938
939 }
940
941 //__________________________________________________________________
942 void AliAnaGammaDirect::Print(const Option_t * opt) const
943 {
944   
945   //Print some relevant parameters set for the analysis
946   if(! opt)
947     return;
948   
949   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
950   
951   printf("Min Gamma pT       =     %2.2f\n",  GetMinPt()) ;
952   printf("Max Gamma pT       =     %3.2f\n",  GetMaxPt()) ;
953
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) ;
960   
961   if(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]) ;
966     printf("    \n") ;
967     
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]) ;
972     
973     printf("    \n") ;
974     
975     printf(" pT fractions         =    \n") ;
976     for(Int_t i = 0; i < fNPtThresFrac; i++)
977       printf("   %2.2f;",  fPtFractions[i]) ;
978     
979   }  
980   
981   printf("    \n") ;
982   
983
984
985 //____________________________________________________________________________
986 Bool_t  AliAnaGammaDirect::SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom){
987    //Select cluster depending on its pid and acceptance selections
988    
989    //Skip matched clusters with tracks
990    if(calo->GetNTracksMatched() > 0) return kFALSE ;
991    
992    //Check PID
993    calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
994    Int_t pdg = AliCaloPID::kPhoton;   
995    if(IsCaloPIDOn()){
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
1000      else
1001        pdg = GetCaloPID()->GetPdg(fDetector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
1002      
1003      if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
1004      
1005      //If it does not pass pid, skip
1006      if(pdg != AliCaloPID::kPhoton) return kFALSE ;
1007    }
1008    
1009    //Check acceptance selection
1010    if(IsFidutialCutOn()){
1011      Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
1012      if(! in ) return kFALSE ;
1013    }
1014    
1015    if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1016    
1017    return kTRUE;
1018    
1019  }