change all printf's to AliDebug/AliInfo/AliWarning
[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---:") ;
101   parList+=onePar ;     
102   snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)",
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   //Initialize the parameters of the analysis.
442   AddToHistogramsName("AnaPhotonConvInCalo_");
443   
444   fMassCut                = 0.03; //30 MeV
445   fRemoveConvertedPair    = kFALSE;
446   fAddConvertedPairsToAOD = kFALSE;
447         
448 }
449
450 //_________________________________________________
451 void  AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()
452 {
453   //Do conversion photon analysis and fill aods
454   
455   //Loop on stored AOD photons
456   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
457   AliDebug(1,Form("AOD branch entries %d", naod));
458   
459   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
460   Bool_t * indexConverted = new Bool_t[naod];
461   for (Int_t i = 0; i < naod; i++) indexConverted[i] = kFALSE;
462   
463   for(Int_t iaod = 0; iaod < naod ; iaod++)
464   {
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       {
488         calo->SetTagged(kFALSE);
489         id2 = calo2->GetCaloLabel(0);
490         Float_t asymmetry = TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E());
491         Float_t dPhi      = (calo->Momentum())->Phi()-(calo2->Momentum())->Phi();
492         Float_t dEta      = (calo->Momentum())->Eta()-(calo2->Momentum())->Eta();
493         
494         //...............................................
495         //Fill few histograms with kinematics of the pair
496         //FIXME, move all this to MakeAnalysisFillHistograms ...
497         
498         fhConvDeltaEta   ->Fill( pairM, dPhi      );
499         fhConvDeltaPhi   ->Fill( pairM, dEta      );
500         fhConvAsym       ->Fill( pairM, asymmetry );
501         fhConvDeltaEtaPhi->Fill( dEta , dPhi      );
502         fhConvPt         ->Fill( pairM, (calo->Momentum())->Pt()+(calo2->Momentum())->Pt());
503         
504         //Estimate conversion distance, T. Awes, M. Ivanov
505         //Under the assumption that the pair has zero mass, and that each electron
506         //of the pair has the same momentum, they will each have the same bend radius
507         //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
508         //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
509         //R = E/1.5 [cm].  Under these assumptions, the distance from the conversion
510         //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
511         //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
512         //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
513         //the clusters.
514         
515         TObjArray * clusters    = 0;
516         if(calo->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
517         else                                 clusters = GetPHOSClusters();
518         
519         
520         Int_t iclus = -1;
521         AliVCluster *cluster1 = FindCluster(clusters,calo ->GetCaloLabel(0),iclus);
522         AliVCluster *cluster2 = FindCluster(clusters,calo2->GetCaloLabel(0),iclus);
523         
524         Float_t pos1[3];
525         cluster1->GetPosition(pos1);
526         Float_t pos2[3];
527         cluster2->GetPosition(pos2);
528         Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
529                                         (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
530                                         (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
531         
532         Float_t convDist  = TMath::Sqrt(calo->E() *clustDist*0.01/0.15);
533         Float_t convDist2 = TMath::Sqrt(calo2->E()*clustDist*0.01/0.15);
534         //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,calo->E(),convDist,calo2->E(),convDist2);
535         AliDebug(2,Form("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",
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         {
548           fhConvDistEtaCutEta ->Fill(calo->Eta(), convDist );
549           fhConvDistEtaCutEta ->Fill(calo2->Eta(),convDist2);
550           fhConvDistEnCutEta  ->Fill(calo->E(),   convDist );
551           fhConvDistEnCutEta  ->Fill(calo2->E(),  convDist2);
552           fhConvDistMassCutEta->Fill(pairM, convDist );
553           //mass cut
554           if(pairM<0.01)
555           {//10 MeV
556             fhConvDistEtaCutMass ->Fill(calo ->Eta(), convDist );
557             fhConvDistEtaCutMass ->Fill(calo2->Eta(), convDist2);
558             fhConvDistEnCutMass  ->Fill(calo ->E(),   convDist );
559             fhConvDistEnCutMass  ->Fill(calo2->E(),   convDist2);
560             // asymmetry cut
561             if(asymmetry<0.1)
562             {
563               fhConvDistEtaCutAsy ->Fill(calo ->Eta(), convDist );
564               fhConvDistEtaCutAsy ->Fill(calo2->Eta(), convDist2);
565               fhConvDistEnCutAsy  ->Fill(calo ->E(),   convDist );
566               fhConvDistEnCutAsy  ->Fill(calo2->E(),   convDist2);
567             }//asymmetry cut
568           }//mass cut
569         }//dEta cut
570         
571         //...............................................
572         //Select pairs in a eta-phi window
573         if(TMath::Abs(dEta) < fConvDEtaCut    &&
574            TMath::Abs(dPhi) < fConvDPhiMaxCut &&
575            TMath::Abs(dPhi) > fConvDPhiMinCut &&
576            asymmetry        < fConvAsymCut       )
577         {
578           indexConverted[iaod] = kTRUE;
579           indexConverted[jaod] = kTRUE;
580           bConverted           = kTRUE;
581         }
582         //printf("Accepted? %d\n",bConverted);
583         //...........................................
584         //Fill more histograms, simulated data
585         //FIXME, move all this to MakeAnalysisFillHistograms ...
586         if(IsDataMC())
587         {
588           //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
589           Int_t ancPDG    = 0;
590           Int_t ancStatus = 0;
591           Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(cluster1->GetLabel(), cluster2->GetLabel(),
592                                                                       GetReader(), ancPDG, ancStatus, fMomentum, fProdVertex);
593           
594           // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
595           //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
596           
597           Int_t tag1 = calo ->GetTag();
598           Int_t tag2 = calo2->GetTag();
599           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion))
600           {
601             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1)
602             {
603               fhConvDeltaEtaMCConversion   ->Fill( pairM, dEta      );
604               fhConvDeltaPhiMCConversion   ->Fill( pairM, dPhi      );
605               fhConvAsymMCConversion       ->Fill( pairM, asymmetry );
606               fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi      );
607               fhConvPtMCConversion         ->Fill( pairM, calo->Pt()+calo2->Pt());
608               fhConvDispersionMCConversion ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
609               fhConvM02MCConversion        ->Fill( cluster1->GetM02(), cluster2->GetM02());
610               fhConvDistMCConversion       ->Fill( convDist , fProdVertex.Mag() );
611               fhConvDistMCConversion       ->Fill( convDist2, fProdVertex.Mag() );
612               
613               if(dEta<0.05 && pairM<0.01 && asymmetry<0.1)
614               {
615                 fhConvDistMCConversionCuts->Fill( convDist , fProdVertex.Mag() );
616                 fhConvDistMCConversionCuts->Fill( convDist2, fProdVertex.Mag() );
617               }
618               
619             }
620           }
621           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiNeutron))
622           {
623             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1)
624             {
625               fhConvDeltaEtaMCAntiNeutron    ->Fill( pairM, dEta      );
626               fhConvDeltaPhiMCAntiNeutron    ->Fill( pairM, dPhi      );
627               fhConvAsymMCAntiNeutron        ->Fill( pairM, asymmetry );
628               fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi      );
629               fhConvPtMCAntiNeutron          ->Fill( pairM, calo->Pt()+calo2->Pt());
630               fhConvDispersionMCAntiNeutron  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
631               fhConvM02MCAntiNeutron         ->Fill( cluster1->GetM02(), cluster2->GetM02());
632             }
633           }
634           else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiProton))
635           {
636             if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1)
637             {
638               fhConvDeltaEtaMCAntiProton    ->Fill( pairM, dEta      );
639               fhConvDeltaPhiMCAntiProton    ->Fill( pairM, dPhi      );
640               fhConvAsymMCAntiProton        ->Fill( pairM, asymmetry );
641               fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi      );
642               fhConvPtMCAntiProton          ->Fill( pairM, calo->Pt()+calo2->Pt());
643               fhConvDispersionMCAntiProton  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
644               fhConvM02MCAntiProton         ->Fill( cluster1->GetM02(), cluster2->GetM02());
645             }
646           }
647           
648           //Pairs coming from fragmenting pairs.
649           if( ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) )
650           {
651             fhConvDeltaEtaMCString    ->Fill( pairM, dPhi);
652             fhConvDeltaPhiMCString    ->Fill( pairM, dPhi);
653             fhConvAsymMCString        ->Fill( pairM, TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E()) );
654             fhConvDeltaEtaPhiMCString ->Fill( dPhi,  dPhi);
655             fhConvPtMCString          ->Fill( pairM, calo->Pt()+calo2->Pt());
656             fhConvDispersionMCString  ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
657             fhConvM02MCString         ->Fill( cluster1->GetM02(), cluster2->GetM02());
658           }
659           
660         }// Data MC
661         
662         break;
663       }
664       
665     }//Mass loop
666     
667     //..........................................................................................................
668     //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
669     if(bConverted)
670     {
671       //Add to AOD
672       if(fAddConvertedPairsToAOD)
673       {
674         //Create AOD of pair analysis
675         fMomentum = *(calo->Momentum())+*(calo2->Momentum());
676         AliAODPWG4Particle aodpair = AliAODPWG4Particle(fMomentum);
677         aodpair.SetLabel(calo->GetLabel());
678         
679         //printf("Index %d, Id %d\n",iaod, calo->GetID());
680         //Set the indeces of the original caloclusters
681         aodpair.SetCaloLabel(calo->GetCaloLabel(0),id2);
682         aodpair.SetDetectorTag(calo->GetDetectorTag());
683         aodpair.SetIdentifiedParticleType(calo->GetIdentifiedParticleType());
684         aodpair.SetTag(calo ->GetTag());
685         aodpair.SetTagged(kTRUE);
686         //Add AOD with pair object to aod branch
687         AddAODParticle(aodpair);
688         //printf("\t \t both added pair\n");
689       }
690       
691       //Do not add the current calocluster
692       if(!fRemoveConvertedPair)
693       {
694         //printf("TAGGED\n");
695         //Tag this cluster as likely conversion
696         calo->SetTagged(kTRUE);
697       }
698     }//converted pair
699     
700   }// main loop
701   
702   // Remove entries identified as conversion electrons
703   // Revise if this is OK
704   if(fRemoveConvertedPair || fAddConvertedPairsToAOD)
705   {
706     for(Int_t iaod = 0; iaod < naod ; iaod++)
707     {
708       if(indexConverted[iaod])GetOutputAODBranch()->RemoveAt(iaod);
709     }
710     GetOutputAODBranch()->Compress();
711   }
712   
713   delete [] indexConverted;
714   
715   AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
716   
717 }
718
719 //________________________________________________________
720 void  AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms()
721 {
722   //Fill histograms
723   
724   //-------------------------------------------------------------------
725   // Access MC information in stack if requested, check that it exists.
726   AliStack         * stack       = 0x0;
727   TParticle        * primary     = 0x0;
728   TClonesArray     * mcparticles = 0x0;
729   AliAODMCParticle * aodprimary  = 0x0;
730   
731   if(IsDataMC())
732   {
733     if(GetReader()->ReadStack())
734     {
735       stack =  GetMCStack() ;
736       if(!stack)
737         AliFatal("Stack not available, is the MC handler called? STOP");
738     }
739     else if(GetReader()->ReadAODMCParticles())
740     {
741       //Get the list of MC particles
742       mcparticles = GetReader()->GetAODMCParticles();
743       if(!mcparticles)
744         AliFatal("Standard MCParticles not available!");
745     }
746   }// is data and MC
747   
748   //----------------------------------
749   //Loop on stored AOD photons
750   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
751   AliDebug(1,Form("AOD branch entries %d", naod));
752   
753   for(Int_t iaod = 0; iaod < naod ; iaod++)
754   {
755     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
756     //Int_t pdg = ph->GetIdentifiedParticleType();
757     
758     if(ph->IsTagged())
759     {
760       AliDebug(2,Form("ID Photon: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta())) ;
761       //................................
762       //Fill photon histograms
763       Float_t ptcluster  = ph->Pt();
764       Float_t phicluster = ph->Phi();
765       Float_t etacluster = ph->Eta();
766       Float_t ecluster   = ph->E();
767       
768       fhPtPhotonConv->Fill(ptcluster);
769       if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
770       else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
771       
772       //.......................................
773       //Play with the MC data if available
774       if(IsDataMC())
775       {
776         //....................................................................
777         // Access MC information in stack if requested, check that it exists.
778         Int_t label =ph->GetLabel();
779         if(label < 0)
780         {
781           AliDebug(1,Form("*** bad label ***:  label %d", label));
782           continue;
783         }
784         
785         //Float_t eprim   = 0;
786         //Float_t ptprim  = 0;
787         if(GetReader()->ReadStack())
788         {
789           if(label >=  stack->GetNtrack())
790           {
791             AliDebug(1,Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
792             continue ;
793           }
794           
795           primary = stack->Particle(label);
796           if(!primary)
797           {
798             AliDebug(1,Form("*** no primary ***:  label %d", label));
799             continue;
800           }
801           //eprim   = primary->Energy();
802           //ptprim  = primary->Pt();
803           
804         }
805         else if(GetReader()->ReadAODMCParticles())
806         {
807           //Check which is the input
808           if(ph->GetInputFileIndex() == 0)
809           {
810             if(!mcparticles) continue;
811             
812             if(label >=  mcparticles->GetEntriesFast())
813             {
814               AliDebug(2,Form("*** large label ***:  label %d, n tracks %d",label, mcparticles->GetEntriesFast()));
815               continue ;
816             }
817             
818             //Get the particle
819             aodprimary = (AliAODMCParticle*) mcparticles->At(label);
820             
821           }
822           
823           if(!aodprimary)
824           {
825             AliDebug(2,Form("*** no primary ***:  label %d", label));
826             continue;
827           }
828           
829           //eprim   = aodprimary->E();
830           //ptprim  = aodprimary->Pt();
831           
832         }
833         
834         Int_t tag =ph->GetTag();
835         
836         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
837         {
838           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
839           {
840             fhPtConversionTagged ->Fill(ptcluster);
841           }
842         }
843         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
844         {
845           fhPtAntiNeutronTagged ->Fill(ptcluster);
846         }
847         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
848         {
849           fhPtAntiProtonTagged ->Fill(ptcluster);
850         }
851         else
852         {
853           fhPtUnknownTagged ->Fill(ptcluster);
854         }
855         
856       }//Histograms with MC
857     }// tagged by conversion
858   }// aod loop
859   
860 }
861
862
863 //________________________________________________________
864 void AliAnaPhotonConvInCalo::Print(const Option_t * opt) const
865 {
866   //Print some relevant parameters set for the analysis
867   
868   if(! opt)
869     return;
870   
871   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
872   AliAnaCaloTrackCorrBaseClass::Print(" ");
873   
874   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
875   printf("Conversion pair mass cut             = %f\n",fMassCut);
876   printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
877          fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
878   
879   printf("    \n") ;
880         
881