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