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