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