bb78fd1d46816a2260ed30c1a3367f36cd89e13f
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPhotonConvInCalo.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 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
16 //_________________________________________________________________________
17 //
18 // Conversions pairs analysis
19 // Check if cluster comes from a conversion in the material in front of the calorimeter
20 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
21 // Input are selected clusters with AliAnaPhoton
22 //
23 //
24 //-- Author: Gustavo Conesa (LPSC-IN2P3-CNRS)
25 //////////////////////////////////////////////////////////////////////////////
26
27
28 // --- ROOT system --- 
29 #include <TH2F.h>
30 #include <TH3D.h>
31 #include <TClonesArray.h>
32 #include <TObjString.h>
33 #include "TParticle.h"
34 #include "TDatabasePDG.h"
35
36 // --- Analysis system --- 
37 #include "AliAnaPhotonConvInCalo.h" 
38 #include "AliCaloTrackReader.h"
39 #include "AliStack.h"
40 #include "AliCaloPID.h"
41 #include "AliMCAnalysisUtils.h"
42 #include "AliFiducialCut.h"
43 #include "AliVCluster.h"
44 #include "AliAODMCParticle.h"
45
46 ClassImp(AliAnaPhotonConvInCalo)
47
48 //________________________________________
49 AliAnaPhotonConvInCalo::AliAnaPhotonConvInCalo() : 
50 AliAnaCaloTrackCorrBaseClass(),   
51 fRemoveConvertedPair(kFALSE), 
52 fAddConvertedPairsToAOD(kFALSE), 
53 fMassCut(0),                  
54 fConvAsymCut(1.),                  fConvDEtaCut(2.),
55 fConvDPhiMinCut(-1.),              fConvDPhiMaxCut(7.), 
56 fMomentum(),                       fProdVertex(),
57 // Histograms
58 fhPtPhotonConv(0),                 fhEtaPhiPhotonConv(0),          fhEtaPhi05PhotonConv(0),
59 fhConvDeltaEta(0),                 fhConvDeltaPhi(0),              fhConvDeltaEtaPhi(0), 
60 fhConvAsym(0),                     fhConvPt(0),
61 fhConvDistEta(0),                  fhConvDistEn(0),                fhConvDistMass(0),     
62 fhConvDistEtaCutEta(0),            fhConvDistEnCutEta(0),          fhConvDistMassCutEta(0),
63 fhConvDistEtaCutMass(0),           fhConvDistEnCutMass(0), 
64 fhConvDistEtaCutAsy(0),            fhConvDistEnCutAsy(0),
65
66 // MC histograms
67 fhPtConversionTagged(0),           fhPtAntiNeutronTagged(0),       
68 fhPtAntiProtonTagged(0),           fhPtUnknownTagged(0),
69
70 fhConvDeltaEtaMCConversion(0),     fhConvDeltaPhiMCConversion(0),  fhConvDeltaEtaPhiMCConversion(0),
71 fhConvAsymMCConversion(0),         fhConvPtMCConversion(0),           
72 fhConvDispersionMCConversion(0),   fhConvM02MCConversion(0),
73
74 fhConvDeltaEtaMCAntiNeutron(0),    fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0), 
75 fhConvAsymMCAntiNeutron(0),        fhConvPtMCAntiNeutron(0), 
76 fhConvDispersionMCAntiNeutron(0),  fhConvM02MCAntiNeutron(0),
77 fhConvDeltaEtaMCAntiProton(0),     fhConvDeltaPhiMCAntiProton(0),  fhConvDeltaEtaPhiMCAntiProton(0),  
78 fhConvAsymMCAntiProton(0),         fhConvPtMCAntiProton(0),  
79 fhConvDispersionMCAntiProton(0),   fhConvM02MCAntiProton(0),
80 fhConvDeltaEtaMCString(0),         fhConvDeltaPhiMCString(0),      fhConvDeltaEtaPhiMCString(0),      
81 fhConvAsymMCString(0),             fhConvPtMCString(0),      
82 fhConvDispersionMCString(0),       fhConvM02MCString(0),
83 fhConvDistMCConversion(0),         fhConvDistMCConversionCuts(0)
84 {
85   //default ctor
86   
87   //Initialize parameters
88   InitParameters();
89   
90 }
91
92 //_________________________________________________
93 TObjString *  AliAnaPhotonConvInCalo::GetAnalysisCuts()
94 {       
95   //Save parameters used for analysis
96   TString parList ; //this will be list of parameters used for this analysis.
97   const Int_t buffersize = 255;
98   char onePar[buffersize] ;
99   
100   snprintf(onePar,buffersize,"--- AliAnaPhotonConvInCalo---\n") ;
101   parList+=onePar ;     
102   snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
103            fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
104   parList+=onePar ; 
105   
106   return new TObjString(parList) ;
107 }
108
109 //___________________________________________________
110 TList *  AliAnaPhotonConvInCalo::GetCreateOutputObjects()
111 {  
112   // Create histograms to be saved in output file and 
113   // store them in outputContainer
114   TList * outputContainer = new TList() ; 
115   outputContainer->SetName("PhotonConvInCaloHistos") ; 
116         
117   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin(); 
118   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin(); 
119   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();   
120   
121   fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax); 
122   fhPtPhotonConv->SetYTitle("N");
123   fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
124   outputContainer->Add(fhPtPhotonConv) ; 
125   
126   fhEtaPhiPhotonConv  = new TH2F
127   ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
128   fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
129   fhEtaPhiPhotonConv->SetXTitle("#eta");
130   outputContainer->Add(fhEtaPhiPhotonConv) ;
131   if(GetMinPt() < 0.5){
132     fhEtaPhi05PhotonConv  = new TH2F
133     ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax); 
134     fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
135     fhEtaPhi05PhotonConv->SetXTitle("#eta");
136     outputContainer->Add(fhEtaPhi05PhotonConv) ;
137   }
138   
139   fhConvDeltaEta  = new TH2F
140   ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5); 
141   fhConvDeltaEta->SetYTitle("#Delta #eta");
142   fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
143   outputContainer->Add(fhConvDeltaEta) ;
144   
145   fhConvDeltaPhi  = new TH2F
146   ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5); 
147   fhConvDeltaPhi->SetYTitle("#Delta #phi");
148   fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
149   outputContainer->Add(fhConvDeltaPhi) ;
150   
151   fhConvDeltaEtaPhi  = new TH2F
152   ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
153   fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
154   fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
155   outputContainer->Add(fhConvDeltaEtaPhi) ;
156   
157   fhConvAsym  = new TH2F
158   ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1); 
159   fhConvAsym->SetYTitle("Asymmetry");
160   fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
161   outputContainer->Add(fhConvAsym) ;  
162   
163   fhConvPt  = new TH2F
164   ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.); 
165   fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
166   fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
167   outputContainer->Add(fhConvPt) ;
168   
169   fhConvDistEta  = new TH2F
170   ("hConvDistEta","distance to conversion vertex",100,-0.7,0.7,100,0.,5.); 
171   fhConvDistEta->SetXTitle("#eta");
172   fhConvDistEta->SetYTitle(" distance (m)");
173   outputContainer->Add(fhConvDistEta) ;
174   
175   fhConvDistEn  = new TH2F
176   ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,100,0.,5.); 
177   fhConvDistEn->SetXTitle("E (GeV)");
178   fhConvDistEn->SetYTitle(" distance (m)");
179   outputContainer->Add(fhConvDistEn) ;
180   
181   fhConvDistMass  = new TH2F
182   ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,100,0.,5.); 
183   fhConvDistMass->SetXTitle("m (GeV/c^2)");
184   fhConvDistMass->SetYTitle(" distance (m)");
185   outputContainer->Add(fhConvDistMass) ;
186   
187   fhConvDistEtaCutEta  = new TH2F
188   ("hConvDistEtaCutEta","distance to conversion vertex, dEta < 0.05",100,-0.7,0.7,100,0.,5.); 
189   fhConvDistEtaCutEta->SetXTitle("#eta");
190   fhConvDistEtaCutEta->SetYTitle(" distance (m)");
191   outputContainer->Add(fhConvDistEtaCutEta) ;
192   
193   fhConvDistEnCutEta  = new TH2F
194   ("hConvDistEnCutEta","distance to conversion vertex, dEta < 0.05",nptbins,ptmin,ptmax,100,0.,5.); 
195   fhConvDistEnCutEta->SetXTitle("E (GeV)");
196   fhConvDistEnCutEta->SetYTitle(" distance (m)");
197   outputContainer->Add(fhConvDistEnCutEta) ;
198   
199   fhConvDistMassCutEta  = new TH2F
200   ("hConvDistMassCutEta","distance to conversion vertex, dEta < 0.05",100,0,fMassCut,100,0.,5.); 
201   fhConvDistMassCutEta->SetXTitle("m (GeV/c^2)");
202   fhConvDistMassCutEta->SetYTitle(" distance (m)");
203   outputContainer->Add(fhConvDistMassCutEta) ;
204   
205   fhConvDistEtaCutMass  = new TH2F
206   ("hConvDistEtaCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",100,-0.7,0.7,100,0.,5.); 
207   fhConvDistEtaCutMass->SetXTitle("#eta");
208   fhConvDistEtaCutMass->SetYTitle(" distance (m)");
209   outputContainer->Add(fhConvDistEtaCutMass) ;
210   
211   fhConvDistEnCutMass  = new TH2F
212   ("hConvDistEnCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",nptbins,ptmin,ptmax,100,0.,5.); 
213   fhConvDistEnCutMass->SetXTitle("E (GeV)");
214   fhConvDistEnCutMass->SetYTitle(" distance (m)");
215   outputContainer->Add(fhConvDistEnCutMass) ;
216   
217   fhConvDistEtaCutAsy  = new TH2F
218   ("hConvDistEtaCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",100,-0.7,0.7,100,0.,5.); 
219   fhConvDistEtaCutAsy->SetXTitle("#eta");
220   fhConvDistEtaCutAsy->SetYTitle(" distance (m)");
221   outputContainer->Add(fhConvDistEtaCutAsy) ;
222   
223   fhConvDistEnCutAsy  = new TH2F
224   ("hConvDistEnCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",nptbins,ptmin,ptmax,100,0.,5.); 
225   fhConvDistEnCutAsy->SetXTitle("E (GeV)");
226   fhConvDistEnCutAsy->SetYTitle(" distance (m)");
227   outputContainer->Add(fhConvDistEnCutAsy) ;
228   
229   if(IsDataMC()){
230     
231     fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
232     fhPtConversionTagged->SetYTitle("N");
233     fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
234     outputContainer->Add(fhPtConversionTagged) ; 
235     
236     
237     fhPtAntiNeutronTagged  = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
238     fhPtAntiNeutronTagged->SetYTitle("N");
239     fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
240     outputContainer->Add(fhPtAntiNeutronTagged) ; 
241     
242     fhPtAntiProtonTagged  = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
243     fhPtAntiProtonTagged->SetYTitle("N");
244     fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
245     outputContainer->Add(fhPtAntiProtonTagged) ; 
246     
247     fhPtUnknownTagged  = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
248     fhPtUnknownTagged->SetYTitle("N");
249     fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
250     outputContainer->Add(fhPtUnknownTagged) ;     
251     
252     fhConvDeltaEtaMCConversion  = new TH2F
253     ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5); 
254     fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
255     fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
256     outputContainer->Add(fhConvDeltaEtaMCConversion) ;
257     
258     fhConvDeltaPhiMCConversion  = new TH2F
259     ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5); 
260     fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
261     fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
262     outputContainer->Add(fhConvDeltaPhiMCConversion) ;
263     
264     fhConvDeltaEtaPhiMCConversion  = new TH2F
265     ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
266     fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
267     fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
268     outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
269     
270     fhConvAsymMCConversion  = new TH2F
271     ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1); 
272     fhConvAsymMCConversion->SetYTitle("Asymmetry");
273     fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
274     outputContainer->Add(fhConvAsymMCConversion) ;
275     
276     fhConvPtMCConversion  = new TH2F
277     ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.); 
278     fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
279     fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
280     outputContainer->Add(fhConvPtMCConversion) ;    
281     
282     fhConvDispersionMCConversion  = new TH2F
283     ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.); 
284     fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
285     fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
286     outputContainer->Add(fhConvDispersionMCConversion) ;   
287     
288     fhConvM02MCConversion  = new TH2F
289     ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
290     fhConvM02MCConversion->SetYTitle("M02 cluster 1");
291     fhConvM02MCConversion->SetXTitle("M02 cluster 2");
292     outputContainer->Add(fhConvM02MCConversion) ;           
293     
294     fhConvDeltaEtaMCAntiNeutron  = new TH2F
295     ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5); 
296     fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
297     fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
298     outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
299     
300     fhConvDeltaPhiMCAntiNeutron  = new TH2F
301     ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5); 
302     fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
303     fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
304     outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
305     
306     fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
307     ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
308     fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
309     fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
310     outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;    
311     
312     fhConvAsymMCAntiNeutron  = new TH2F
313     ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1); 
314     fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
315     fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
316     outputContainer->Add(fhConvAsymMCAntiNeutron) ;
317     
318     fhConvPtMCAntiNeutron  = new TH2F
319     ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.); 
320     fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
321     fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
322     outputContainer->Add(fhConvPtMCAntiNeutron) ;    
323     
324     fhConvDispersionMCAntiNeutron  = new TH2F
325     ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.); 
326     fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
327     fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
328     outputContainer->Add(fhConvDispersionMCAntiNeutron) ;       
329     
330     fhConvM02MCAntiNeutron  = new TH2F
331     ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
332     fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
333     fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
334     outputContainer->Add(fhConvM02MCAntiNeutron) ;  
335     
336     fhConvDeltaEtaMCAntiProton  = new TH2F
337     ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5); 
338     fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
339     fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
340     outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
341     
342     fhConvDeltaPhiMCAntiProton  = new TH2F
343     ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5); 
344     fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
345     fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
346     outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
347     
348     fhConvDeltaEtaPhiMCAntiProton  = new TH2F
349     ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
350     fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
351     fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
352     outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;    
353     
354     fhConvAsymMCAntiProton  = new TH2F
355     ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1); 
356     fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
357     fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
358     outputContainer->Add(fhConvAsymMCAntiProton) ;
359     
360     fhConvPtMCAntiProton  = new TH2F
361     ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.); 
362     fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
363     fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
364     outputContainer->Add(fhConvPtMCAntiProton) ;
365     
366     fhConvDispersionMCAntiProton  = new TH2F
367     ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.); 
368     fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
369     fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
370     outputContainer->Add(fhConvDispersionMCAntiProton) ;       
371     
372     fhConvM02MCAntiProton  = new TH2F
373     ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
374     fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
375     fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
376     outputContainer->Add(fhConvM02MCAntiProton) ;       
377     
378     fhConvDeltaEtaMCString  = new TH2F
379     ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5); 
380     fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
381     fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
382     outputContainer->Add(fhConvDeltaEtaMCString) ;
383     
384     fhConvDeltaPhiMCString  = new TH2F
385     ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5); 
386     fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
387     fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
388     outputContainer->Add(fhConvDeltaPhiMCString) ;
389     
390     fhConvDeltaEtaPhiMCString  = new TH2F
391     ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
392     fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
393     fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
394     outputContainer->Add(fhConvDeltaEtaPhiMCString) ;    
395     
396     fhConvAsymMCString  = new TH2F
397     ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1); 
398     fhConvAsymMCString->SetYTitle("Asymmetry");
399     fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
400     outputContainer->Add(fhConvAsymMCString) ;
401     
402     fhConvPtMCString  = new TH2F
403     ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.); 
404     fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
405     fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
406     outputContainer->Add(fhConvPtMCString) ;
407     
408     fhConvDispersionMCString  = new TH2F
409     ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
410     fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
411     fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
412     outputContainer->Add(fhConvDispersionMCString) ;       
413     
414     fhConvM02MCString  = new TH2F
415     ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
416     fhConvM02MCString->SetYTitle("M02 cluster 1");
417     fhConvM02MCString->SetXTitle("M02 cluster 2");
418     outputContainer->Add(fhConvM02MCString) ; 
419     
420     fhConvDistMCConversion  = new TH2F
421     ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",100,0.,5.,100,0.,5.); 
422     fhConvDistMCConversion->SetYTitle("distance");
423     fhConvDistMCConversion->SetXTitle("vertex R");
424     outputContainer->Add(fhConvDistMCConversion) ; 
425     
426     fhConvDistMCConversionCuts  = new TH2F
427     ("hConvDistMCConversionCuts","calculated conversion distance vs real vertes for MC conversion, deta < 0.05, m < 10 MeV, asym < 0.1",100,0.,5.,100,0.,5.); 
428     fhConvDistMCConversionCuts->SetYTitle("distance");
429     fhConvDistMCConversionCuts->SetXTitle("vertex R");
430     outputContainer->Add(fhConvDistMCConversionCuts) ; 
431         
432   }
433   
434   return outputContainer ;
435
436 }
437
438 //_______________________________________
439 void AliAnaPhotonConvInCalo::InitParameters()
440 {
441   
442   //Initialize the parameters of the analysis.
443   AddToHistogramsName("AnaPhotonConvInCalo_");
444   
445   fMassCut                = 0.03; //30 MeV
446   fRemoveConvertedPair    = kFALSE;
447   fAddConvertedPairsToAOD = kFALSE;
448         
449 }
450
451 //_____________________________________________
452 void  AliAnaPhotonConvInCalo::MakeAnalysisFillAOD() 
453 {
454   //Do conversion photon analysis and fill aods
455   
456   //Loop on stored AOD photons
457   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
458   if(GetDebug() > 0) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillAOD() - aod branch entries %d\n", naod);
459   
460   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
461   Bool_t * indexConverted = new Bool_t[naod];
462   for (Int_t i = 0; i < naod; i++) indexConverted[i] = kFALSE;
463         
464   for(Int_t iaod = 0; iaod < naod ; iaod++){
465     AliAODPWG4Particle* calo =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
466     
467     Bool_t bConverted = kFALSE;
468     Int_t id2 = -1;
469     
470     //Check if set previously as converted couple, if so skip its use.
471     if (indexConverted[iaod]) continue;
472     
473     // Second cluster loop
474     AliAODPWG4Particle* calo2 = 0;
475     for(Int_t jaod = iaod + 1 ; jaod < naod ; jaod++)
476     {
477       //Check if set previously as converted couple, if so skip its use.
478       if (indexConverted[jaod]) continue;
479       //printf("Check Conversion indeces %d and %d\n",iaod,jaod);
480       calo2 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(jaod));
481       
482       //................................................
483       //Get mass of pair, if small, take this pair.
484       Float_t pairM     = calo->GetPairMass(calo2);
485       //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
486       if(pairM < fMassCut){  
487         calo->SetTagged(kFALSE);
488         id2 = calo2->GetCaloLabel(0);
489         Float_t asymmetry = TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E());
490         Float_t dPhi      = (calo->Momentum())->Phi()-(calo2->Momentum())->Phi();
491         Float_t dEta      = (calo->Momentum())->Eta()-(calo2->Momentum())->Eta();  
492         
493         //...............................................
494         //Fill few histograms with kinematics of the pair
495         //FIXME, move all this to MakeAnalysisFillHistograms ...
496         
497         fhConvDeltaEta   ->Fill( pairM, dPhi      );
498         fhConvDeltaPhi   ->Fill( pairM, dEta      );
499         fhConvAsym       ->Fill( pairM, asymmetry );
500         fhConvDeltaEtaPhi->Fill( dEta , dPhi      );
501         fhConvPt         ->Fill( pairM, (calo->Momentum())->Pt()+(calo2->Momentum())->Pt());          
502         
503         //Estimate conversion distance, T. Awes, M. Ivanov
504         //Under the assumption that the pair has zero mass, and that each electron 
505         //of the pair has the same momentum, they will each have the same bend radius 
506         //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m]. 
507         //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,  
508         //R = E/1.5 [cm].  Under these assumptions, the distance from the conversion 
509         //point to the MCEal can be related to the separation distance, L=2y, on the MCEal 
510         //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as 
511         //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between 
512         //the clusters.
513         
514         TObjArray * clusters    = 0; 
515         if(calo->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
516         else                                 clusters = GetPHOSClusters();
517         
518         
519         Int_t iclus = -1;
520         AliVCluster *cluster1 = FindCluster(clusters,calo ->GetCaloLabel(0),iclus); 
521         AliVCluster *cluster2 = FindCluster(clusters,calo2->GetCaloLabel(0),iclus); 
522
523         Float_t pos1[3];
524         cluster1->GetPosition(pos1); 
525         Float_t pos2[3];
526         cluster2->GetPosition(pos2); 
527         Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
528                                         (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
529                                         (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
530         
531         Float_t convDist  = TMath::Sqrt(calo->E() *clustDist*0.01/0.15);
532         Float_t convDist2 = TMath::Sqrt(calo2->E()*clustDist*0.01/0.15);
533         //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,calo->E(),convDist,calo2->E(),convDist2);
534         if(GetDebug() > 2)
535           printf("AliAnaPhotonConvInCalo::MakeAnalysisFillAOD(): Pair with mass %2.3f < %2.3f, %1.2f < dPhi %2.2f < %2.2f, dEta %f < %2.2f, asymmetry %2.2f< %2.2f; \n    cluster1 id %d, e %2.3f  SM %d, eta %2.3f, phi %2.3f ; \n    cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f\n",
536                  pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
537                  calo->GetCaloLabel(0),calo->E(),GetCaloUtils()->GetModuleNumber(calo,GetReader()->GetInputEvent()), calo->Eta(), calo->Phi(),
538                  id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2,GetReader()->GetInputEvent()),calo2->Eta(), calo2->Phi());
539         
540         fhConvDistEta ->Fill(calo ->Eta(),convDist );
541         fhConvDistEta ->Fill(calo2->Eta(),convDist2);
542         fhConvDistEn  ->Fill(calo ->E(),  convDist );
543         fhConvDistEn  ->Fill(calo2->E(),  convDist2);        
544         fhConvDistMass->Fill(pairM, convDist );
545         //dEta cut
546         if(dEta<0.05){
547           fhConvDistEtaCutEta ->Fill(calo->Eta(), convDist );
548           fhConvDistEtaCutEta ->Fill(calo2->Eta(),convDist2);
549           fhConvDistEnCutEta  ->Fill(calo->E(),   convDist );
550           fhConvDistEnCutEta  ->Fill(calo2->E(),  convDist2);        
551           fhConvDistMassCutEta->Fill(pairM, convDist );
552           //mass cut
553           if(pairM<0.01){//10 MeV
554             fhConvDistEtaCutMass ->Fill(calo ->Eta(), convDist );
555             fhConvDistEtaCutMass ->Fill(calo2->Eta(), convDist2);
556             fhConvDistEnCutMass  ->Fill(calo ->E(),   convDist );
557             fhConvDistEnCutMass  ->Fill(calo2->E(),   convDist2);        
558             // asymmetry cut
559             if(asymmetry<0.1){
560               fhConvDistEtaCutAsy ->Fill(calo ->Eta(), convDist );
561               fhConvDistEtaCutAsy ->Fill(calo2->Eta(), convDist2);
562               fhConvDistEnCutAsy  ->Fill(calo ->E(),   convDist );
563               fhConvDistEnCutAsy  ->Fill(calo2->E(),   convDist2); 
564             }//asymmetry cut
565           }//mass cut            
566         }//dEta cut
567         
568         //...............................................
569         //Select pairs in a eta-phi window
570         if(TMath::Abs(dEta) < fConvDEtaCut    && 
571            TMath::Abs(dPhi) < fConvDPhiMaxCut &&
572            TMath::Abs(dPhi) > fConvDPhiMinCut && 
573            asymmetry        < fConvAsymCut       ){
574           indexConverted[iaod] = kTRUE;
575           indexConverted[jaod] = kTRUE; 
576           bConverted           = kTRUE;          
577         }
578         //printf("Accepted? %d\n",bConverted);
579         //...........................................
580         //Fill more histograms, simulated data
581         //FIXME, move all this to MakeAnalysisFillHistograms ...
582         if(IsDataMC()){
583           
584           //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
585           Int_t ancPDG    = 0;
586           Int_t ancStatus = 0;
587           Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(cluster1->GetLabel(), cluster2->GetLabel(), 
588                                                                       GetReader(), ancPDG, ancStatus, fMomentum, fProdVertex);
589           
590           // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
591           //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
592           
593           Int_t tag1 = calo ->GetTag();
594           Int_t tag2 = calo2->GetTag();
595           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion)){
596             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
597               fhConvDeltaEtaMCConversion   ->Fill( pairM, dEta      );
598               fhConvDeltaPhiMCConversion   ->Fill( pairM, dPhi      );
599               fhConvAsymMCConversion       ->Fill( pairM, asymmetry );
600               fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi      );
601               fhConvPtMCConversion         ->Fill( pairM, calo->Pt()+calo2->Pt());
602               fhConvDispersionMCConversion ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
603               fhConvM02MCConversion        ->Fill( cluster1->GetM02(), cluster2->GetM02());
604               fhConvDistMCConversion       ->Fill( convDist , fProdVertex.Mag() );
605               fhConvDistMCConversion       ->Fill( convDist2, fProdVertex.Mag() );
606               
607               if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
608                 fhConvDistMCConversionCuts->Fill( convDist , fProdVertex.Mag() );
609                 fhConvDistMCConversionCuts->Fill( convDist2, fProdVertex.Mag() );
610               }
611               
612             }              
613           }
614           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiNeutron)){
615             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
616               fhConvDeltaEtaMCAntiNeutron    ->Fill( pairM, dEta      );
617               fhConvDeltaPhiMCAntiNeutron    ->Fill( pairM, dPhi      );
618               fhConvAsymMCAntiNeutron        ->Fill( pairM, asymmetry );
619               fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi      );
620               fhConvPtMCAntiNeutron          ->Fill( pairM, calo->Pt()+calo2->Pt());
621               fhConvDispersionMCAntiNeutron  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
622               fhConvM02MCAntiNeutron         ->Fill( cluster1->GetM02(), cluster2->GetM02());
623             }
624           }
625           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiProton)){
626             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
627               fhConvDeltaEtaMCAntiProton    ->Fill( pairM, dEta      );
628               fhConvDeltaPhiMCAntiProton    ->Fill( pairM, dPhi      );
629               fhConvAsymMCAntiProton        ->Fill( pairM, asymmetry );
630               fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi      );
631               fhConvPtMCAntiProton          ->Fill( pairM, calo->Pt()+calo2->Pt());
632               fhConvDispersionMCAntiProton  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
633               fhConvM02MCAntiProton         ->Fill( cluster1->GetM02(), cluster2->GetM02());
634             }
635           }
636           
637           //Pairs coming from fragmenting pairs.
638           if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
639             fhConvDeltaEtaMCString    ->Fill( pairM, dPhi);
640             fhConvDeltaPhiMCString    ->Fill( pairM, dPhi);
641             fhConvAsymMCString        ->Fill( pairM, TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E()) );
642             fhConvDeltaEtaPhiMCString ->Fill( dPhi,  dPhi);
643             fhConvPtMCString          ->Fill( pairM, calo->Pt()+calo2->Pt());
644             fhConvDispersionMCString  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
645             fhConvM02MCString         ->Fill( cluster1->GetM02(), cluster2->GetM02());
646           }
647           
648         }// Data MC
649         
650         break;
651       }
652       
653     }//Mass loop
654     
655     //..........................................................................................................
656     //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
657     if(bConverted)
658     {
659       //Add to AOD
660       if(fAddConvertedPairsToAOD){
661         //Create AOD of pair analysis
662         fMomentum = *(calo->Momentum())+*(calo2->Momentum());
663         AliAODPWG4Particle aodpair = AliAODPWG4Particle(fMomentum);
664         aodpair.SetLabel(calo->GetLabel());
665         
666         //printf("Index %d, Id %d\n",iaod, calo->GetID());
667         //Set the indeces of the original caloclusters  
668         aodpair.SetCaloLabel(calo->GetCaloLabel(0),id2);
669         aodpair.SetDetectorTag(calo->GetDetectorTag());
670         aodpair.SetIdentifiedParticleType(calo->GetIdentifiedParticleType());
671         aodpair.SetTag(calo ->GetTag());
672         aodpair.SetTagged(kTRUE);
673         //Add AOD with pair object to aod branch
674         AddAODParticle(aodpair);
675         //printf("\t \t both added pair\n");
676       }
677       
678       //Do not add the current calocluster
679       if(!fRemoveConvertedPair) 
680       {
681         //printf("TAGGED\n");
682         //Tag this cluster as likely conversion
683         calo->SetTagged(kTRUE);
684       }
685     }//converted pair
686     
687   }// main loop
688   
689   // Remove entries identified as conversion electrons
690   // Revise if this is OK
691   if(fRemoveConvertedPair || fAddConvertedPairsToAOD){
692     for(Int_t iaod = 0; iaod < naod ; iaod++)
693       if(indexConverted[iaod])GetOutputAODBranch()->RemoveAt(iaod);
694     GetOutputAODBranch()->Compress();
695   }
696   
697   delete [] indexConverted;
698         
699   if(GetDebug() > 1) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());  
700   
701 }
702
703 //____________________________________________________
704 void  AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() 
705 {
706   //Fill histograms
707   
708   //-------------------------------------------------------------------
709   // Access MC information in stack if requested, check that it exists. 
710   AliStack         * stack       = 0x0;
711   TParticle        * primary     = 0x0;   
712   TClonesArray     * mcparticles = 0x0;
713   AliAODMCParticle * aodprimary  = 0x0; 
714   
715   if(IsDataMC()){
716     
717     if(GetReader()->ReadStack()){
718       stack =  GetMCStack() ;
719       if(!stack) {
720         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
721         abort();
722       }
723       
724     }
725     else if(GetReader()->ReadAODMCParticles()){
726       
727       //Get the list of MC particles
728       mcparticles = GetReader()->GetAODMCParticles();
729       if(!mcparticles && GetDebug() > 0)
730       {
731         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
732       } 
733     }
734   }// is data and MC
735  
736   //----------------------------------
737   //Loop on stored AOD photons
738   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
739   if(GetDebug() > 0) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
740   
741   for(Int_t iaod = 0; iaod < naod ; iaod++){
742     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
743     //Int_t pdg = ph->GetIdentifiedParticleType();
744     
745     if(ph->IsTagged()){
746       
747       if(GetDebug() > 2) 
748         printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
749       //................................
750       //Fill photon histograms 
751       Float_t ptcluster  = ph->Pt();
752       Float_t phicluster = ph->Phi();
753       Float_t etacluster = ph->Eta();
754       Float_t ecluster   = ph->E();
755       
756       fhPtPhotonConv->Fill(ptcluster);
757       if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
758       else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
759       
760       
761       //.......................................
762       //Play with the MC data if available
763       if(IsDataMC()){
764         
765         
766         //....................................................................
767         // Access MC information in stack if requested, check that it exists.
768         Int_t label =ph->GetLabel();
769         if(label < 0) {
770           if(GetDebug() > 1) printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
771           continue;
772         }
773         
774         //Float_t eprim   = 0;
775         //Float_t ptprim  = 0;
776         if(GetReader()->ReadStack()){
777           
778           if(label >=  stack->GetNtrack()) {
779             if(GetDebug() > 2)  printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
780             continue ;
781           }
782           
783           primary = stack->Particle(label);
784           if(!primary){
785             printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
786             continue;
787           }
788           //eprim   = primary->Energy();
789           //ptprim  = primary->Pt();
790           
791         }
792         else if(GetReader()->ReadAODMCParticles()){
793           //Check which is the input
794           if(ph->GetInputFileIndex() == 0){
795             if(!mcparticles) continue;
796             if(label >=  mcparticles->GetEntriesFast()) {
797               if(GetDebug() > 2)  printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", 
798                                          label, mcparticles->GetEntriesFast());
799               continue ;
800             }
801             //Get the particle
802             aodprimary = (AliAODMCParticle*) mcparticles->At(label);
803             
804           }
805           
806           if(!aodprimary){
807             printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
808             continue;
809           }
810           
811           //eprim   = aodprimary->E();
812           //ptprim  = aodprimary->Pt();
813           
814         }
815         
816         Int_t tag =ph->GetTag();
817         
818         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
819         {
820           
821           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
822           {
823             
824             fhPtConversionTagged ->Fill(ptcluster);
825             
826           }                     
827         }
828         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
829         {
830           
831           fhPtAntiNeutronTagged ->Fill(ptcluster);
832           
833         }
834         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
835         {
836           fhPtAntiProtonTagged ->Fill(ptcluster);
837           
838         } 
839         
840         else {
841           fhPtUnknownTagged ->Fill(ptcluster);
842           
843         }
844         
845       }//Histograms with MC
846     }// tagged by conversion
847   }// aod loop
848   
849 }
850
851
852 //________________________________________________________
853 void AliAnaPhotonConvInCalo::Print(const Option_t * opt) const
854 {
855   //Print some relevant parameters set for the analysis
856   
857   if(! opt)
858     return;
859   
860   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
861   AliAnaCaloTrackCorrBaseClass::Print(" ");
862   
863   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
864   printf("Conversion pair mass cut             = %f\n",fMassCut);
865   printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
866          fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
867   
868   printf("    \n") ;
869         
870