AliPhysics  b76e98e (b76e98e)
AliAnaPhotonConvInCalo.cxx
Go to the documentation of this file.
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 // --- ROOT system ---
17 #include <TH2F.h>
18 #include <TH3D.h>
19 #include <TClonesArray.h>
20 #include <TObjString.h>
21 #include "TDatabasePDG.h"
22 
23 // --- Analysis system ---
24 #include "AliAnaPhotonConvInCalo.h"
25 #include "AliCaloTrackReader.h"
26 #include "AliMCEvent.h"
27 #include "AliCaloPID.h"
28 #include "AliMCAnalysisUtils.h"
29 #include "AliFiducialCut.h"
30 #include "AliVCluster.h"
31 #include "AliVParticle.h"
32 
34 ClassImp(AliAnaPhotonConvInCalo) ;
36 
37 //________________________________________________
39 //________________________________________________
42 fRemoveConvertedPair(kFALSE),
43 fAddConvertedPairsToAOD(kFALSE),
44 fFillClusterConvDistHisto(kFALSE),
45 fMassCut(0), fMassCutTight(0),
46 fConvAsymCut(1.), fConvDEtaCut(2.),
47 fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
48 fMomentum(), fProdVertex(),
49 // Histograms
50 fhPtPhotonConv(0),
51 fhConvDeltaEta(0), fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0),
52 fhConvAsym(0), fhConvPt(0),
53 fhConvDistEta(0), fhConvDistPhi(0), fhConvDistEn(0), fhConvDistMass(0),
54 fhConvDistEtaCutEta(0), fhConvDistPhiCutEta(0), fhConvDistEnCutEta(0), fhConvDistMassCutEta(0),
55 fhConvDistEtaCutMass(0), fhConvDistPhiCutMass(0), fhConvDistEnCutMass(0),
56 fhConvDistEtaCutAsy(0), fhConvDistPhiCutAsy(0), fhConvDistEnCutAsy(0), fhConvDistMassCutAsy(0),
57 fhConvDistEtaCutAll(0), fhConvDistPhiCutAll(0), fhConvDistEnCutAll(0),
58 
59 // MC histograms
60 fhPtConversionTagged(0), fhPtAntiNeutronTagged(0),
61 fhPtAntiProtonTagged(0), fhPtUnknownTagged(0),
62 
63 fhConvDeltaEtaMCConversion(0), fhConvDeltaPhiMCConversion(0), fhConvDeltaEtaPhiMCConversion(0),
64 fhConvAsymMCConversion(0), fhConvPtMCConversion(0),
65 //fhConvDispersionMCConversion(0),
66 fhConvM02MCConversion(0),
67 
68 fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0),
69 fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0),
70 //fhConvDispersionMCAntiNeutron(0),
71 fhConvM02MCAntiNeutron(0),
72 fhConvDeltaEtaMCAntiProton(0), fhConvDeltaPhiMCAntiProton(0), fhConvDeltaEtaPhiMCAntiProton(0),
73 fhConvAsymMCAntiProton(0), fhConvPtMCAntiProton(0),
74 //fhConvDispersionMCAntiProton(0),
75 fhConvM02MCAntiProton(0),
76 fhConvDeltaEtaMCString(0), fhConvDeltaPhiMCString(0), fhConvDeltaEtaPhiMCString(0),
77 fhConvAsymMCString(0), fhConvPtMCString(0),
78 //fhConvDispersionMCString(0),
79 fhConvM02MCString(0),
80 fhConvDistMCConversion(0), fhConvDistMCConversionCuts(0)
81 {
83 
84  for(Int_t ibin = 0; ibin < 6; ibin++)
85  {
86  fhEtaPhiPhotonConv [ibin] = 0;
87  fhEtaPhiPhotonConvPaired[ibin] = 0;
88  fhConvPtMCConversionRcut[ibin] = 0;
89  //fhConvPtRcut [ibin] = 0;
90  }
91 
92 }
93 
94 //_____________________________________________________
96 //_____________________________________________________
98 {
99  TString parList ; //this will be list of parameters used for this analysis.
100  const Int_t buffersize = 255;
101  char onePar[buffersize] ;
102 
103  snprintf(onePar,buffersize,"--- AliAnaPhotonConvInCalo---:") ;
104  parList+=onePar ;
105  snprintf(onePar,buffersize,"Conversion pair mass cut = %1.2f and tight mass cut = %1.2f; ",fMassCut,fMassCutTight);
106  parList+=onePar ;
107  snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f); ",
109  parList+=onePar ;
110  snprintf(onePar,buffersize,"Add conversion pair to AOD = %d, remove partners? %d. ",fAddConvertedPairsToAOD, fRemoveConvertedPair);
111  parList+=onePar ;
112 
113  return new TObjString(parList) ;
114 }
115 
116 //_______________________________________________________
119 //_______________________________________________________
121 {
122  TList * outputContainer = new TList() ;
123  outputContainer->SetName("PhotonConvInCaloHistos") ;
124 
128 
129  Int_t ndist = 250;
130  Int_t mindist = 0;
131  Int_t maxdist = 500;
132 
133  TString region[] = {"ITS","TPC","TRD","TOF","Top EMCal","In EMCal"};
134 
135  fhPtPhotonConv = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
136  fhPtPhotonConv->SetYTitle("N");
137  fhPtPhotonConv->SetXTitle("#it{p}_{T #gamma}(GeV/#it{c})");
138  outputContainer->Add(fhPtPhotonConv) ;
139 
140  Float_t ebin[] = {0.3,0.4,0.5,0.6,0.75,1,2};
141 
142  for(Int_t iebin = 0; iebin < 6; iebin++)
143  {
144  fhEtaPhiPhotonConv[iebin] = new TH2F
145  (Form("hEtaPhiPhotonConv_ebin%d",iebin),Form("pair #eta vs #phi, %2.2f<#it{E}<%2.2f GeV/#it{c}",ebin[iebin],ebin[iebin+1]),
146  netabins,etamin,etamax,nphibins,phimin,phimax);
147  fhEtaPhiPhotonConv[iebin]->SetYTitle("#phi (rad)");
148  fhEtaPhiPhotonConv[iebin]->SetXTitle("#eta");
149  outputContainer->Add(fhEtaPhiPhotonConv[iebin]) ;
150 
151  fhEtaPhiPhotonConvPaired[iebin] = new TH2F
152  (Form("hEtaPhiPhotonConvPaired_ebin%d",iebin),Form("cluster #eta vs #phi, %2.2f<#it{E}<%2.2f GeV/#it{c}",ebin[iebin],ebin[iebin+1]),
153  netabins,etamin,etamax,nphibins,phimin,phimax);
154  fhEtaPhiPhotonConvPaired[iebin]->SetYTitle("#phi (rad)");
155  fhEtaPhiPhotonConvPaired[iebin]->SetXTitle("#eta");
156  outputContainer->Add(fhEtaPhiPhotonConvPaired[iebin]) ;
157  }
158 
159  fhConvDeltaEta = new TH2F
160  ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
161  fhConvDeltaEta->SetYTitle("#Delta #eta");
162  fhConvDeltaEta->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
163  outputContainer->Add(fhConvDeltaEta) ;
164 
165  fhConvDeltaPhi = new TH2F
166  ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
167  fhConvDeltaPhi->SetYTitle("#Delta #phi");
168  fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
169  outputContainer->Add(fhConvDeltaPhi) ;
170 
171  fhConvDeltaEtaPhi = new TH2F
172  ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
173  fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
174  fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
175  outputContainer->Add(fhConvDeltaEtaPhi) ;
176 
177  fhConvAsym = new TH2F
178  ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
179  fhConvAsym->SetYTitle("Asymmetry");
180  fhConvAsym->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
181  outputContainer->Add(fhConvAsym) ;
182 
183  fhConvPt = new TH2F
184  ("hConvPt","#it{p}_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
185  fhConvPt->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
186  fhConvPt->SetXTitle("Pair #it{M} (GeV/#it{c}^{2})");
187  outputContainer->Add(fhConvPt) ;
188 
189 // for(Int_t iR = 0; iR < 6; iR++)
190 // {
191 // fhConvPtRcut[iR] = new TH2F
192 // (Form("hConvPt_R%d",iR),
193 // Form("#it{p}_{T} of selected conversion pairs in %s?",region[iR].Data()),
194 // 100,0,fMassCut,100,0.,10.);
195 // fhConvPtRcut[iR]->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
196 // fhConvPtRcut[iR]->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
197 // outputContainer->Add(fhConvPtRcut[iR]) ;
198 // }
199 
201  {
202  fhConvDistEta = new TH2F
203  ("hConvDistEta","distance to conversion vertex",netabins,etamin,etamax,ndist,mindist,maxdist);
204  fhConvDistEta->SetXTitle("#eta");
205  fhConvDistEta->SetYTitle(" distance (cm)");
206  outputContainer->Add(fhConvDistEta) ;
207 
208  fhConvDistPhi = new TH2F
209  ("hConvDistPhi","distance to conversion vertex",nphibins,phimin,phimax,ndist,mindist,maxdist);
210  fhConvDistPhi->SetXTitle("#phi (rad)");
211  fhConvDistPhi->SetYTitle(" distance (cm)");
212  outputContainer->Add(fhConvDistPhi) ;
213 
214  fhConvDistEn = new TH2F
215  ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,ndist,mindist,maxdist);
216  fhConvDistEn->SetXTitle("#it{E} (GeV)");
217  fhConvDistEn->SetYTitle(" distance (cm)");
218  outputContainer->Add(fhConvDistEn) ;
219 
220  fhConvDistMass = new TH2F
221  ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,ndist,mindist,maxdist);
222  fhConvDistMass->SetXTitle("#it{M} (GeV/#it{c}^{2})");
223  fhConvDistMass->SetYTitle(" distance (cm)");
224  outputContainer->Add(fhConvDistMass) ;
225 
227  ("hConvDistEtaCutEta",Form("distance to conversion vertex, #Delta #eta < %2.2f",fConvDEtaCut),netabins,etamin,etamax,ndist,mindist,maxdist);
228  fhConvDistEtaCutEta->SetXTitle("#eta");
229  fhConvDistEtaCutEta->SetYTitle(" distance (cm)");
230  outputContainer->Add(fhConvDistEtaCutEta) ;
231 
233  ("hConvDistPhiCutEta",Form("distance to conversion vertex, #Delta #eta < %2.2f",fConvDEtaCut),nphibins,phimin,phimax,ndist,mindist,maxdist);
234  fhConvDistPhiCutEta->SetXTitle("#phi (rad)");
235  fhConvDistPhiCutEta->SetYTitle(" distance (cm)");
236  outputContainer->Add(fhConvDistPhiCutEta) ;
237 
239  ("hConvDistEnCutEta",Form("distance to conversion vertex, #Delta #eta < %2.2f",fConvDEtaCut),nptbins,ptmin,ptmax,ndist,mindist,maxdist);
240  fhConvDistEnCutEta->SetXTitle("#it{E} (GeV)");
241  fhConvDistEnCutEta->SetYTitle(" distance (cm)");
242  outputContainer->Add(fhConvDistEnCutEta) ;
243 
245  ("hConvDistMassCutEta",Form("distance to conversion vertex, #Delta #eta < %2.2f",fConvDEtaCut),100,0,fMassCut,ndist,mindist,maxdist);
246  fhConvDistMassCutEta->SetXTitle("#it{M} (GeV/#it{c}^{2})");
247  fhConvDistMassCutEta->SetYTitle(" distance (cm)");
248  outputContainer->Add(fhConvDistMassCutEta) ;
249 
251  ("hConvDistEtaCutMass",Form("distance to conversion vertex, #it{M} < %0.2f MeV/#it{c}^{2}",fMassCutTight),netabins,etamin,etamax,ndist,mindist,maxdist);
252  fhConvDistEtaCutMass->SetXTitle("#eta");
253  fhConvDistEtaCutMass->SetYTitle(" distance (cm)");
254  outputContainer->Add(fhConvDistEtaCutMass) ;
255 
257  ("hConvDistPhiCutMass",Form("distance to conversion vertex, #it{M} < %0.2f MeV/#it{c}^{2}",fMassCutTight),nphibins,phimin,phimax,ndist,mindist,maxdist);
258  fhConvDistPhiCutMass->SetXTitle("#phi (rad)");
259  fhConvDistPhiCutMass->SetYTitle(" distance (cm)");
260  outputContainer->Add(fhConvDistPhiCutMass) ;
261 
263  ("hConvDistEnCutMass",Form("distance to conversion vertex, #it{M} < %0.2f MeV/#it{c}^{2}",fMassCutTight),nptbins,ptmin,ptmax,ndist,mindist,maxdist);
264  fhConvDistEnCutMass->SetXTitle("#it{E} (GeV)");
265  fhConvDistEnCutMass->SetYTitle(" distance (cm)");
266  outputContainer->Add(fhConvDistEnCutMass) ;
267 
269  ("hConvDistEtaCutAsy",Form("distance to conversion vertex, #it{A} < %2.2f",fConvAsymCut),netabins,etamin,etamax,ndist,mindist,maxdist);
270  fhConvDistEtaCutAsy->SetXTitle("#eta");
271  fhConvDistEtaCutAsy->SetYTitle(" distance (cm)");
272  outputContainer->Add(fhConvDistEtaCutAsy) ;
273 
275  ("hConvDistPhiCutAsy",Form("distance to conversion vertex, #it{A} < %2.2f",fConvAsymCut),nphibins,phimin,phimax,ndist,mindist,maxdist);
276  fhConvDistPhiCutAsy->SetXTitle("#phi (rad)");
277  fhConvDistPhiCutAsy->SetYTitle(" distance (cm)");
278  outputContainer->Add(fhConvDistPhiCutAsy) ;
279 
281  ("hConvDistEnCutAsy",Form("distance to conversion vertex, #it{A} < %2.2f",fConvAsymCut),nptbins,ptmin,ptmax,ndist,mindist,maxdist);
282  fhConvDistEnCutAsy->SetXTitle("#it{E} (GeV)");
283  fhConvDistEnCutAsy->SetYTitle(" distance (cm)");
284  outputContainer->Add(fhConvDistEnCutAsy) ;
285 
287  ("hConvDistMassCutAsy",Form("distance to conversion vertex, #it{A} < %2.2f",fConvDEtaCut),100,0,fMassCut,ndist,mindist,maxdist);
288  fhConvDistMassCutAsy->SetXTitle("#it{M} (GeV/#it{c}^{2})");
289  fhConvDistMassCutEta->SetYTitle(" distance (cm)");
290  outputContainer->Add(fhConvDistMassCutAsy) ;
291 
293  ("hConvDistEtaCutAll",Form("distance to conversion vertex, #it{M} < %0.2f MeV/#it{c}^{2}, #Delta #eta < %2.2f, #it{A} < %2.2f",
294  fConvDEtaCut,fMassCutTight,fConvAsymCut),netabins,etamin,etamax,ndist,mindist,maxdist);
295  fhConvDistEtaCutAll->SetXTitle("#eta");
296  fhConvDistEtaCutAll->SetYTitle(" distance (cm)");
297  outputContainer->Add(fhConvDistEtaCutAll) ;
298 
300  ("hConvDistPhiCutAll",Form("distance to conversion vertex, #it{M} < %0.2f MeV/#it{c}^{2}, #Delta #eta < %2.2f, #it{A} < %2.2f",
301  fConvDEtaCut,fMassCutTight,fConvAsymCut),nphibins,phimin,phimax,ndist,mindist,maxdist);
302  fhConvDistPhiCutAll->SetXTitle("#phi (rad)");
303  fhConvDistPhiCutAll->SetYTitle(" distance (cm)");
304  outputContainer->Add(fhConvDistPhiCutAll) ;
305 
307  ("hConvDistEnCutAll",Form("distance to conversion vertex, #it{M} < %0.2f MeV/#it{c}^{2}, #Delta #eta < %2.2f, #it{A} < %2.2f",
308  fConvDEtaCut,fMassCutTight,fConvAsymCut),nptbins,ptmin,ptmax,ndist,mindist,maxdist);
309  fhConvDistEnCutAll->SetXTitle("#it{E} (GeV)");
310  fhConvDistEnCutAll->SetYTitle(" distance (cm)");
311  outputContainer->Add(fhConvDistEnCutAll) ;
312  }
313 
314  if(IsDataMC())
315  {
316  fhPtConversionTagged = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
317  fhPtConversionTagged->SetYTitle("N");
318  fhPtConversionTagged->SetXTitle("#it{p}_{T #gamma}(GeV/#it{c})");
319  outputContainer->Add(fhPtConversionTagged) ;
320 
321  fhPtAntiNeutronTagged = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
322  fhPtAntiNeutronTagged->SetYTitle("N");
323  fhPtAntiNeutronTagged->SetXTitle("#it{p}_{T #gamma}(GeV/#it{c})");
324  outputContainer->Add(fhPtAntiNeutronTagged) ;
325 
326  fhPtAntiProtonTagged = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
327  fhPtAntiProtonTagged->SetYTitle("N");
328  fhPtAntiProtonTagged->SetXTitle("#it{p}_{T #gamma}(GeV/#it{c})");
329  outputContainer->Add(fhPtAntiProtonTagged) ;
330 
331  fhPtUnknownTagged = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
332  fhPtUnknownTagged->SetYTitle("N");
333  fhPtUnknownTagged->SetXTitle("#it{p}_{T #gamma}(GeV/#it{c})");
334  outputContainer->Add(fhPtUnknownTagged) ;
335 
337  {
339  ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",ndist,mindist,maxdist,ndist,mindist,maxdist);
340  fhConvDistMCConversion->SetXTitle("distance");
341  fhConvDistMCConversion->SetYTitle("vertex R");
342  outputContainer->Add(fhConvDistMCConversion) ;
343 
345  ("hConvDistMCConversionCuts",
346  Form("calculated conversion distance vs real vertes for MC conversion, #Delta #eta < %2.2f, #it{M} < 10 MeV/#it{c}^{2}, #it{A} < %2.2f",
348  ndist,mindist,maxdist,ndist,mindist,maxdist);
349  fhConvDistMCConversionCuts->SetXTitle("distance");
350  fhConvDistMCConversionCuts->SetYTitle("vertex R");
351  outputContainer->Add(fhConvDistMCConversionCuts) ;
352  }
353 
355  ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
356  fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
357  fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
358  outputContainer->Add(fhConvDeltaEtaMCConversion) ;
359 
361  ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
362  fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
363  fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
364  outputContainer->Add(fhConvDeltaPhiMCConversion) ;
365 
367  ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
368  fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
369  fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
370  outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
371 
373  ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
374  fhConvAsymMCConversion->SetYTitle("Asymmetry");
375  fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
376  outputContainer->Add(fhConvAsymMCConversion) ;
377 
379  ("hConvPtMCConversion","#it{p}_{T} of selected conversion pairs from real conversions",100,0,fMassCut,nptbins,ptmin,ptmax);
380  fhConvPtMCConversion->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
381  fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
382  outputContainer->Add(fhConvPtMCConversion) ;
383 
384  for(Int_t iR = 0; iR < 6; iR++)
385  {
387  (Form("hConvPtMCConversion_R%d",iR),
388  Form("#it{p}_{T} of selected conversion pairs from real conversions in %s",region[iR].Data()),
389  100,0,fMassCut,nptbins,ptmin,ptmax);
390  fhConvPtMCConversionRcut[iR]->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
391  fhConvPtMCConversionRcut[iR]->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
392  outputContainer->Add(fhConvPtMCConversionRcut[iR]) ;
393  }
394 
395 // fhConvDispersionMCConversion = new TH2F
396 // ("hConvDispersionMCConversion","#it{p}_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
397 // fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
398 // fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
399 // outputContainer->Add(fhConvDispersionMCConversion) ;
400 
402  ("hConvM02MCConversion","#it{p}_{T} of selected conversion pairs from real conversion",100,0.,1.,100,0.,1.);
403  fhConvM02MCConversion->SetYTitle("M02 cluster 1");
404  fhConvM02MCConversion->SetXTitle("M02 cluster 2");
405  outputContainer->Add(fhConvM02MCConversion) ;
406 
408  ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
409  fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
410  fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
411  outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
412 
414  ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
415  fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
416  fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
417  outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
418 
420  ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
421  fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
422  fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
423  outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
424 
426  ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
427  fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
428  fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
429  outputContainer->Add(fhConvAsymMCAntiNeutron) ;
430 
432  ("hConvPtMCAntiNeutron","#it{p}_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,nptbins,ptmin,ptmax);
433  fhConvPtMCAntiNeutron->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
434  fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
435  outputContainer->Add(fhConvPtMCAntiNeutron) ;
436 
437 // fhConvDispersionMCAntiNeutron = new TH2F
438 // ("hConvDispersionMCAntiNeutron","#it{p}_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
439 // fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
440 // fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
441 // outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
442 
444  ("hConvM02MCAntiNeutron","#it{p}_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
445  fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
446  fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
447  outputContainer->Add(fhConvM02MCAntiNeutron) ;
448 
450  ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
451  fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
452  fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
453  outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
454 
456  ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
457  fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
458  fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
459  outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
460 
462  ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
463  fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
464  fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
465  outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
466 
468  ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
469  fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
470  fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
471  outputContainer->Add(fhConvAsymMCAntiProton) ;
472 
474  ("hConvPtMCAntiProton","#it{p}_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,nptbins,ptmin,ptmax);
475  fhConvPtMCAntiProton->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
476  fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
477  outputContainer->Add(fhConvPtMCAntiProton) ;
478 
479 // fhConvDispersionMCAntiProton = new TH2F
480 // ("hConvDispersionMCAntiProton","#it{p}_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
481 // fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
482 // fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
483 // outputContainer->Add(fhConvDispersionMCAntiProton) ;
484 
486  ("hConvM02MCAntiProton","#it{p}_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
487  fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
488  fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
489  outputContainer->Add(fhConvM02MCAntiProton) ;
490 
492  ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
493  fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
494  fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
495  outputContainer->Add(fhConvDeltaEtaMCString) ;
496 
498  ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
499  fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
500  fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
501  outputContainer->Add(fhConvDeltaPhiMCString) ;
502 
504  ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
505  fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
506  fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
507  outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
508 
510  ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
511  fhConvAsymMCString->SetYTitle("Asymmetry");
512  fhConvAsymMCString->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
513  outputContainer->Add(fhConvAsymMCString) ;
514 
515  fhConvPtMCString = new TH2F
516  ("hConvPtMCString","#it{p}_{T} of selected conversion pairs from string",100,0,fMassCut,nptbins,ptmin,ptmax);
517  fhConvPtMCString->SetYTitle("Pair #it{p}_{T} (GeV/#it{c})");
518  fhConvPtMCString->SetXTitle("Pair Mass (GeV/#it{c}^{2})");
519  outputContainer->Add(fhConvPtMCString) ;
520 
521 // fhConvDispersionMCString = new TH2F
522 // ("hConvDispersionMCString","#it{p}_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
523 // fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
524 // fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
525 // outputContainer->Add(fhConvDispersionMCString) ;
526 
527  fhConvM02MCString = new TH2F
528  ("hConvM02MCString","#it{p}_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
529  fhConvM02MCString->SetYTitle("M02 cluster 1");
530  fhConvM02MCString->SetXTitle("M02 cluster 2");
531  outputContainer->Add(fhConvM02MCString) ;
532  }
533 
534  return outputContainer ;
535 }
536 
537 //___________________________________________
539 //___________________________________________
541 {
542  AddToHistogramsName("AnaPhotonConvInCalo_");
543 
544  fConvAsymCut = 0.2 ;
545  fMassCut = 0.04 ; // 40 MeV/c^2
546  fMassCutTight = 0.03 ; // 30 MeV/c^2
547  fConvDEtaCut = 0.04 ;
548  fConvDPhiMinCut = 0 ;
549  fConvDPhiMaxCut = 0.05 ;
550  fRemoveConvertedPair = kFALSE;
551  fAddConvertedPairsToAOD = kFALSE;
552 }
553 
554 //________________________________________________________
557 //________________________________________________________
559 {
560  // Loop on stored AOD photons
561  Int_t naod = GetOutputAODBranch()->GetEntriesFast();
562  AliDebug(1,Form("AOD branch entries %d", naod));
563 
564  // List to be used in conversion analysis, to tag the cluster as candidate for conversion
565  Bool_t * indexConverted = new Bool_t[naod];
566  for (Int_t i = 0; i < naod; i++) indexConverted[i] = kFALSE;
567 
568  for(Int_t iaod = 0; iaod < naod ; iaod++)
569  {
571 
572  //Bool_t bConverted = kFALSE;
573 
574  // Check if set previously as converted couple, if so skip its use.
575  //if (indexConverted[iaod]) continue;
576 
577  // Second cluster loop
578  AliCaloTrackParticle* calo2 = 0;
579  for(Int_t jaod = iaod + 1 ; jaod < naod ; jaod++)
580  {
581  // Check if set previously as converted couple, if so skip its use.
582  //if (indexConverted[jaod]) continue;
583 
584  //printf("Check Conversion indeces %d and %d\n",iaod,jaod);
585  calo2 = (AliCaloTrackParticle*) (GetOutputAODBranch()->At(jaod));
586 
587  //
588  // Add both clusters
589  //
590  fMomentum = *(calo->GetMomentum())+*(calo2->GetMomentum());
591  Float_t ptConvPair = fMomentum.Pt();
592  Float_t phiConvPair = fMomentum.Phi();
593  Float_t etaConvPair = fMomentum.Eta();
594  //Float_t eConvPair = fMomentum.E();
595  Float_t pairM = calo->GetPairMass(calo2);
596  //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
597 
598  //................................................
599  // Get mass of pair, if small, take this pair.
600  if(pairM < fMassCut)
601  {
602  Float_t phi = calo->Phi();
603  if(phi < 0) phi +=TMath::TwoPi();
604  Float_t phi2 = calo2->Phi();
605  if(phi2 < 0) phi2+=TMath::TwoPi();
606 
607  calo->SetTagged(kFALSE);
608  Float_t asymmetry = TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E());
609  Float_t dPhi = phi-phi2;
610  Float_t dEta = calo->Eta() - calo2->Eta();
611 
612  //...............................................
613  // Fill few histograms with kinematics of the pair
614  // FIXME, move all this to MakeAnalysisFillHistograms ...
615 
616  fhConvDeltaEta ->Fill( pairM, dEta , GetEventWeight());
617  fhConvDeltaPhi ->Fill( pairM, dPhi , GetEventWeight());
618  fhConvAsym ->Fill( pairM, asymmetry , GetEventWeight());
619  fhConvDeltaEtaPhi->Fill( dEta , dPhi , GetEventWeight());
620  fhConvPt ->Fill( pairM, ptConvPair, GetEventWeight());
621 
622  Float_t convDist = -1;
623  Float_t convDist2 = -1;
625  {
626  // Estimate conversion distance, T. Awes, M. Ivanov
627  // Under the assumption that the pair has zero mass, and that each electron
628  // of the pair has the same momentum, they will each have the same bend radius
629  // given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
630  // With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
631  // R = E/1.5 [cm]. Under these assumptions, the distance from the conversion
632  // point to the EMCal can be related to the separation distance, L=2y, on the EMCal
633  // as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
634  // d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
635  // the clusters.
636 
637  TObjArray * clusters = 0;
638  if(calo->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
639  else clusters = GetPHOSClusters ();
640 
641  Int_t iclus = -1;
642  AliVCluster *cluster1 = FindCluster(clusters,calo ->GetCaloLabel(0),iclus);
643  AliVCluster *cluster2 = FindCluster(clusters,calo2->GetCaloLabel(0),iclus);
644 
645  Float_t pos1[3];
646  cluster1->GetPosition(pos1);
647  Float_t pos2[3];
648  cluster2->GetPosition(pos2);
649  Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
650  (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
651  (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
652 
653  Float_t convDist = TMath::Sqrt(calo ->E()*clustDist*0.01/0.15) * 100.; // cm
654  Float_t convDist2 = TMath::Sqrt(calo2->E()*clustDist*0.01/0.15) * 100.; // cm
655  //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,calo->E(),convDist,calo2->E(),convDist2);
656  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"
657  " cluster1 id %d, e %2.3f SM %d, eta %2.3f, phi %2.3f ; \n"
658  " cluster2 id %d, e %2.3f, SM %d, eta %2.3f, phi %2.3f \n",
659  pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
660  calo ->GetCaloLabel(0), calo ->E(), GetCaloUtils()->GetModuleNumber(calo ,GetReader()->GetInputEvent()), calo ->Eta(), phi,
661  calo2->GetCaloLabel(0), calo2->E(), GetCaloUtils()->GetModuleNumber(calo2,GetReader()->GetInputEvent()), calo2->Eta(), phi2));
662 
663 
664 // Int_t convDistR1 = -1;
665 // if ( convDist > 430-75. ) convDistR1 = 0;
666 // else if ( convDist < 430-275. ) convDistR1 = 1;
667 // else if ( convDist < 430-375. ) convDistR1 = 2;
668 // else if ( convDist < 430-400. ) convDistR1 = 3;
669 // else if ( convDist < 0 ) convDistR1 = 4;
670 // else convDistR1 = 5;
671 //
672 // Int_t convDistR2 = -1;
673 // if ( convDist2 > 430-75. ) convDistR2 = 0;
674 // else if ( convDist2 < 430-275. ) convDistR2 = 1;
675 // else if ( convDist2 < 430-375. ) convDistR2 = 2;
676 // else if ( convDist2 < 430-400. ) convDistR2 = 3;
677 // else if ( convDist2 < 0 ) convDistR2 = 4;
678 // else convDistR2 = 5;
679 //
680 // if(convDistR1 >= 0)
681 // fhConvPtRcut[convDistR1]->Fill( pairM, ptConvPair, GetEventWeight());
682 // if(convDistR2 >= 0)
683 // fhConvPtRcut[convDistR2]->Fill( pairM, ptConvPair, GetEventWeight());
684 
685  fhConvDistEta ->Fill(calo ->Eta(),convDist , GetEventWeight());
686  fhConvDistEta ->Fill(calo2->Eta(),convDist2, GetEventWeight());
687  fhConvDistPhi ->Fill(phi , convDist , GetEventWeight());
688  fhConvDistPhi ->Fill(phi2, convDist2, GetEventWeight());
689  fhConvDistEn ->Fill(calo ->E(), convDist , GetEventWeight());
690  fhConvDistEn ->Fill(calo2->E(), convDist2, GetEventWeight());
691  fhConvDistMass->Fill(pairM, convDist , GetEventWeight());
692  fhConvDistMass->Fill(pairM, convDist2, GetEventWeight());
693 
694  //dEta cut
695  if(TMath::Abs(dEta) < fConvDEtaCut)
696  {
697  fhConvDistEtaCutEta ->Fill(calo ->Eta(), convDist , GetEventWeight());
698  fhConvDistEtaCutEta ->Fill(calo2->Eta(), convDist2, GetEventWeight());
699  fhConvDistPhiCutEta ->Fill(phi , convDist , GetEventWeight());
700  fhConvDistPhiCutEta ->Fill(phi2, convDist2, GetEventWeight());
701  fhConvDistEnCutEta ->Fill(calo ->E(), convDist , GetEventWeight());
702  fhConvDistEnCutEta ->Fill(calo2->E(), convDist2, GetEventWeight());
703  fhConvDistMassCutEta->Fill(pairM, convDist , GetEventWeight());
704  fhConvDistMassCutEta->Fill(pairM, convDist2, GetEventWeight());
705  }
706 
707  //mass cut
708  if(pairM < fMassCutTight)
709  {
710  fhConvDistEtaCutMass ->Fill(calo ->Eta(), convDist , GetEventWeight());
711  fhConvDistEtaCutMass ->Fill(calo2->Eta(), convDist2, GetEventWeight());
712  fhConvDistPhiCutMass ->Fill(phi , convDist , GetEventWeight());
713  fhConvDistPhiCutMass ->Fill(phi2, convDist2, GetEventWeight());
714  fhConvDistEnCutMass ->Fill(calo ->E(), convDist , GetEventWeight());
715  fhConvDistEnCutMass ->Fill(calo2->E(), convDist2, GetEventWeight());
716  }
717 
718  // asymmetry cut
719  if(asymmetry < fConvAsymCut)
720  {
721  fhConvDistEtaCutAsy ->Fill(calo ->Eta(), convDist , GetEventWeight());
722  fhConvDistEtaCutAsy ->Fill(calo2->Eta(), convDist2, GetEventWeight());
723  fhConvDistPhiCutAsy ->Fill(phi , convDist , GetEventWeight());
724  fhConvDistPhiCutAsy ->Fill(phi2, convDist2, GetEventWeight());
725  fhConvDistEnCutAsy ->Fill(calo ->E(), convDist , GetEventWeight());
726  fhConvDistEnCutAsy ->Fill(calo2->E(), convDist2, GetEventWeight());
727  fhConvDistMassCutAsy->Fill(pairM, convDist , GetEventWeight());
728  fhConvDistMassCutAsy->Fill(pairM, convDist2, GetEventWeight());
729  }// asymmetry cut
730 
731  // All cuts
732  if(TMath::Abs(dEta) < fConvDEtaCut && pairM < fMassCutTight && asymmetry < fConvAsymCut)
733  {
734  fhConvDistEtaCutAll ->Fill(calo ->Eta(), convDist , GetEventWeight());
735  fhConvDistEtaCutAll ->Fill(calo2->Eta(), convDist2, GetEventWeight());
736  fhConvDistPhiCutAll ->Fill(phi , convDist , GetEventWeight());
737  fhConvDistPhiCutAll ->Fill(phi2, convDist2, GetEventWeight());
738  fhConvDistEnCutAll ->Fill(calo ->E(), convDist , GetEventWeight());
739  fhConvDistEnCutAll ->Fill(calo2->E(), convDist2, GetEventWeight());
740  }
741  }
742 
743  //...........................................
744  // Fill more histograms, simulated data
745  Int_t ancPDG = 0;
746  Int_t ancStatus = 0;
747  Int_t ancLabel =-1;
748  if(IsDataMC())
749  {
750  // Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
751 
752  ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
753  GetMC(), ancPDG, ancStatus, fMomentum, fProdVertex);
754 
755  // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
756  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
757 
758  Int_t tag1 = calo ->GetTag();
759  Int_t tag2 = calo2->GetTag();
760 // Float_t disp1 = cluster1->GetDispersion();
761 // Float_t disp2 = cluster2->GetDispersion();
762  Float_t l0cl1 = calo ->GetM02();
763  Float_t l0cl2 = calo2->GetM02();
764 
766  {
767  if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1)
768  {
769  fhConvDeltaEtaMCConversion ->Fill( pairM, dEta , GetEventWeight());
770  fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi , GetEventWeight());
771  fhConvAsymMCConversion ->Fill( pairM, asymmetry , GetEventWeight());
772  fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi , GetEventWeight());
773  fhConvPtMCConversion ->Fill( pairM, ptConvPair, GetEventWeight());
774  //fhConvDispersionMCConversion ->Fill( disp1, disp2 , GetEventWeight());
775  fhConvM02MCConversion ->Fill( l0cl1, l0cl2 , GetEventWeight());
776 
777  Float_t prodR = TMath::Sqrt(fProdVertex.X()*fProdVertex.X()+fProdVertex.Y()*fProdVertex.Y());
778 
780  {
781  fhConvDistMCConversion ->Fill( convDist , prodR, GetEventWeight());
782  fhConvDistMCConversion ->Fill( convDist2, prodR, GetEventWeight());
783  }
784 
785  Int_t convR = -1;
786  if ( prodR < 75. ) convR = 0;
787  else if ( prodR < 275. ) convR = 1;
788  else if ( prodR < 375. ) convR = 2;
789  else if ( prodR < 400. ) convR = 3;
790  else if ( prodR < 430. ) convR = 4;
791  else convR = 5;
792 
793  if ( convR >= 0 )
794  fhConvPtMCConversionRcut[convR]->Fill( pairM, ptConvPair, GetEventWeight());
795 
796  if ( fFillClusterConvDistHisto && dEta < fConvDEtaCut && pairM < fMassCutTight && asymmetry < fConvAsymCut )
797  {
798  fhConvDistMCConversionCuts->Fill( convDist , prodR, GetEventWeight());
799  fhConvDistMCConversionCuts->Fill( convDist2, prodR, GetEventWeight());
800  }
801 
802  }
803  }
804  else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiNeutron))
805  {
806  if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1)
807  {
808  fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta , GetEventWeight());
809  fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi , GetEventWeight());
810  fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry , GetEventWeight());
811  fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi , GetEventWeight());
812  fhConvPtMCAntiNeutron ->Fill( pairM, ptConvPair, GetEventWeight());
813  //fhConvDispersionMCAntiNeutron ->Fill( disp1, disp2 , GetEventWeight());
814  fhConvM02MCAntiNeutron ->Fill( l0cl1, l0cl2 , GetEventWeight());
815  }
816  }
817  else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiProton))
818  {
819  if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1)
820  {
821  fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta , GetEventWeight());
822  fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi , GetEventWeight());
823  fhConvAsymMCAntiProton ->Fill( pairM, asymmetry , GetEventWeight());
824  fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi , GetEventWeight());
825  fhConvPtMCAntiProton ->Fill( pairM, ptConvPair, GetEventWeight());
826  //fhConvDispersionMCAntiProton ->Fill( disp1, disp2 , GetEventWeight());
827  fhConvM02MCAntiProton ->Fill( l0cl1, l0cl2 , GetEventWeight());
828  }
829  }
830 
831  // Pairs coming from fragmenting pairs.
832  if( ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) )
833  {
834  fhConvDeltaEtaMCString ->Fill( pairM, dEta , GetEventWeight());
835  fhConvDeltaPhiMCString ->Fill( pairM, dPhi , GetEventWeight());
836  fhConvAsymMCString ->Fill( pairM, asymmetry , GetEventWeight());
837  fhConvDeltaEtaPhiMCString ->Fill( dEta, dPhi , GetEventWeight());
838  fhConvPtMCString ->Fill( pairM, ptConvPair, GetEventWeight());
839  //fhConvDispersionMCString ->Fill( disp1, disp2 , GetEventWeight());
840  fhConvM02MCString ->Fill( l0cl1, l0cl2 , GetEventWeight());
841  }
842 
843  }// Data MC
844 
845  //...............................................
846  //...............................................
847  // Select pairs in a eta/phi/asymmetry windows
848  //
849  if(TMath::Abs(dEta) < fConvDEtaCut &&
850  TMath::Abs(dPhi) < fConvDPhiMaxCut &&
851  TMath::Abs(dPhi) > fConvDPhiMinCut &&
852  asymmetry < fConvAsymCut )
853  {
854  indexConverted[iaod] = kTRUE;
855  indexConverted[jaod] = kTRUE;
856  //bConverted = kTRUE;
857 
858  Float_t ebin[] = {0.3,0.4,0.5,0.6,0.75,1,2};
859  Int_t bin1 = -1;
860  Int_t bin2 = -1;
861  Int_t bin12 = -1;
862  for(Int_t iebin = 0; iebin < 6; iebin++)
863  {
864  if( calo ->Pt() > ebin[iebin] && calo ->Pt() <= ebin[iebin+1] ) bin1 = iebin ;
865  if( calo2->Pt() > ebin[iebin] && calo2->Pt() <= ebin[iebin+1] ) bin2 = iebin ;
866  if( ptConvPair > ebin[iebin] && ptConvPair <= ebin[iebin+1] ) bin12 = iebin ;
867  }
868 
869  if(bin1 > -1) fhEtaPhiPhotonConvPaired[bin1 ]->Fill(calo ->Eta(), phi , GetEventWeight());
870  if(bin2 > -1) fhEtaPhiPhotonConvPaired[bin2 ]->Fill(calo2->Eta(), phi2 , GetEventWeight());
871  if(bin12 > -1) fhEtaPhiPhotonConv [bin12]->Fill(etaConvPair , phiConvPair, GetEventWeight());
872 
873  fhPtPhotonConv->Fill(ptConvPair, GetEventWeight());
874 
875  if(IsDataMC())
876  {
877  //....................................................................
878  // Access MC information in stack if requested, check that it exists.
879 
880  Int_t label =calo->GetLabel();
881  if ( label < 0 )
882  {
883  AliDebug(1,Form("*** bad label ***: label %d", label));
884  continue;
885  }
886 
887  if ( label >= GetMC()->GetNumberOfTracks() )
888  {
889  AliDebug(1,Form("*** large label ***: label %d, n tracks %d", label, GetMC()->GetNumberOfTracks()));
890  continue ;
891  }
892 
893  //Float_t eprim = 0;
894  //Float_t ptprim = 0;
895 
896  // Get the particle
897  AliVParticle* primary = GetMC()->GetTrack(label);
898 
899  if(!primary)
900  {
901  AliDebug(2,Form("*** no primary ***: label %d", label));
902  continue;
903  }
904  //eprim = aodprimary->E();
905  //ptprim = aodprimary->Pt();
906 
907  Int_t tag = calo ->GetTag();
908  if(ancLabel >=0 )
909  {
910  if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
911  {
913  {
914  fhPtConversionTagged ->Fill(ptConvPair, GetEventWeight());
915  }
916  }
917  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
918  {
919  fhPtAntiNeutronTagged ->Fill(ptConvPair, GetEventWeight());
920  }
921  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
922  {
923  fhPtAntiProtonTagged ->Fill(ptConvPair, GetEventWeight());
924  }
925  }
926  else
927  {
928  fhPtUnknownTagged ->Fill(ptConvPair, GetEventWeight());
929  }
930  } // MC
931 
932  //..........................................................................................................
933  // Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
934  // Add to AOD
936  {
937  // Create AOD of pair analysis
939  aodpair.SetLabel(calo->GetLabel());
940 
941  // Set the indeces of the original caloclusters
942  aodpair.SetCaloLabel(calo->GetCaloLabel(0),calo2->GetCaloLabel(0));
943  aodpair.SetDetectorTag(calo->GetDetectorTag());
945  aodpair.SetTag(calo->GetTag());
946  aodpair.SetTagged(kTRUE);
947 
948  // Add AOD with pair object to aod branch
949  AddAODParticle(aodpair);
950  //printf("added pair: naod %d new %d\n",naod,GetOutputAODBranch()->GetEntriesFast());
951  }
952 
953  // Do not add the current calocluster
955  {
956  //printf("TAGGED\n");
957 
958  // Tag this cluster as likely conversion
959  calo ->SetTagged(kTRUE);
960  calo2->SetTagged(kTRUE);
961  }
962  } // selected converted pair
963  } // mass cut
964 
965  } // Mass loop
966 
967  }// main loop
968 
969  // Remove entries identified as conversion electrons
970  // Revise if this is OK
972  {
973  for(Int_t iaod = 0; iaod < naod ; iaod++)
974  {
975  if(indexConverted[iaod]) GetOutputAODBranch()->RemoveAt(iaod);
976  }
977  GetOutputAODBranch()->Compress();
978  }
979 
980  delete [] indexConverted;
981 
982  AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
983 }
984 
985 
986 //____________________________________________________________
988 //____________________________________________________________
990 {
991  if(! opt)
992  return;
993 
994  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
996 
997  printf("Add conversion pair to AOD = %d, removed? %d\n",fAddConvertedPairsToAOD, fRemoveConvertedPair);
998  printf("Conversion pair mass cut = %1.2f and %1.2f\n",fMassCut,fMassCutTight);
999  printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
1001 
1002  printf(" \n") ;
1003 }
1004 
Float_t GetHistoPtMax() const
TH2F * fhConvDeltaEtaPhi
! Small mass photons, correlation in phi and eta
TH2F * fhConvDeltaEtaMCAntiNeutron
! Small mass cluster pairs, correlation in eta, origin of both clusters is anti neutron ...
virtual Double_t Eta() const
TH2F * fhConvM02MCAntiProton
! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is anti proton ...
TH2F * fhConvDistMassCutEta
! Approx distance to vertex vs Mass, dEta < fConvDEtaCut
virtual Double_t Pt() const
Float_t GetHistoPtMin() const
Float_t fConvDPhiMaxCut
Select conversion pairs when dphi of pair smaller than cut.
Float_t fConvDPhiMinCut
Select conversion pairs when dphi of pair lager than cut.
virtual void AddAODParticle(AliCaloTrackParticle part)
virtual void AddToHistogramsName(TString add)
TH2F * fhConvM02MCAntiNeutron
! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is anti neutron ...
Definition: External.C:236
TH2F * fhConvDistEnCutAll
! Approx distance to vertex vs energy, dEta < fConvDEtacut, M < 20 MeV/c^2, A < fConvAsymCut ...
TH2F * fhConvDistPhi
! Approx distance to vertex vs azimuth
TH2F * fhConvDistMCConversionCuts
! Calculated conversion distance vs real distance to vertex
Bool_t fRemoveConvertedPair
Remove conversion pairs.
Float_t fConvDEtaCut
Select conversion pairs when deta of pair smaller than cut.
TH2F * fhConvM02MCConversion
! Small mass cluster pairs, m02 of cluster 1 vs cluster 2
virtual Int_t GetLabel() const
TH2F * fhConvDeltaEtaMCConversion
! Small mass cluster pairs, correlation in eta, origin of both clusters is conversion ...
TH2F * fhConvPt
! Small mass photons, pT of pair
Double_t GetPairMass(AliCaloTrackParticle *p) const
TH2F * fhConvDeltaPhiMCConversion
! Small mass cluster pairs, correlation in phi, origin of both clusters is conversion ...
TH2F * fhConvDistEta
! Approx distance to vertex vs rapidity
Int_t GetHistoPhiBins() const
virtual void SetDetectorTag(UInt_t d)
TH1F * fhPtAntiNeutronTagged
! Number of identified gamma from AntiNeutrons gamma, tagged as conversion
virtual TClonesArray * GetOutputAODBranch() const
virtual UInt_t GetDetectorTag() const
TH2F * fhConvDeltaEta
! Small mass photons, correlation in eta
virtual void SetIdentifiedParticleType(Int_t pdg)
Float_t GetHistoPhiMin() const
TH2F * fhConvAsymMCString
! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is string ...
virtual Double_t E() const
TH2F * fhConvPtMCConversionRcut[6]
! Small mass cluster pairs, pt of pair, origin of both clusters is conversion, for different producti...
TH2F * fhConvPtMCConversion
! Small mass cluster pairs, pt of pair, origin of both clusters is conversion
TH2F * fhConvAsym
! Small mass photons, correlation in energy asymmetry
Int_t GetModuleNumber(AliCaloTrackParticle *particle, AliVEvent *inputEvent) const
Get the EMCAL/PHOS module number that corresponds to this particle.
const Double_t etamin
TH2F * fhConvDeltaEtaMCString
! Small mass cluster pairs, correlation in eta, origin of both clusters is string ...
TH2F * fhConvDistEtaCutAsy
! Approx distance to vertex vs rapidity, A < fConvAsymCut
Base class for CaloTrackCorr analysis algorithms.
TH2F * fhConvDeltaPhi
! Small mass photons, correlation in phi
TH2F * fhConvDeltaPhiMCAntiProton
! Small mass cluster pairs, correlation in phi, origin of both clusters is anti proton ...
TH1F * fhPtUnknownTagged
! Number of identified gamma from unknown, tagged as conversion
TH2F * fhConvDistEtaCutEta
! Approx distance to vertex vs rapidity, dEta < fConvDEtaCut
void InitParameters()
Initialize the parameters of the analysis.
TObjString * GetAnalysisCuts()
Save parameters used for analysis.
TH1F * fhPtPhotonConv
! Number of identified photon vs transverse momentum
TH2F * fhConvDistEnCutAsy
! Approx distance to vertex vs energy, A < fConvAsymCut
int Int_t
Definition: External.C:63
Container for input particle information on CaloTrackCorr package.
TH2F * fhConvPtMCAntiNeutron
! Small mass cluster pairs, pt of pair, origin of both clusters is anti neutron
virtual AliHistogramRanges * GetHistogramRanges()
TH2F * fhConvDistEnCutEta
! Approx distance to vertex vs Energy, dEta < fConvDEtaCut
float Float_t
Definition: External.C:68
TH2F * fhConvDeltaEtaPhiMCAntiProton
! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is anti proton ...
TH2F * fhConvPtMCString
! Small mass cluster pairs, pt of pairs, origin of both clusters is string
const Double_t ptmax
TH2F * fhConvDistMass
! Approx distance to vertex vs Mass
TH2F * fhEtaPhiPhotonConv[6]
! Pseudorapidity vs Phi of identified photon conv pair for 6 transverse momentum bins ...
TH2F * fhConvDeltaPhiMCAntiNeutron
! Small mass cluster pairs, correlation in phi, origin of both clusters is anti neutron ...
TH1F * fhPtConversionTagged
! Number of identified gamma from Conversion , tagged as conversion
TH2F * fhConvDistPhiCutMass
! Approx distance to vertex vs azimuth, M < 20 MeV/c^2
TLorentzVector * GetMomentum() const
virtual AliCalorimeterUtils * GetCaloUtils() const
TH2F * fhConvDeltaEtaPhiMCConversion
! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is conversion ...
AliAnaPhotonConvInCalo()
Default constructor. Initialize parameters.
TH2F * fhEtaPhiPhotonConvPaired[6]
! Pseudorapidity vs Phi of identified photon conv leg for 6 transverse momentum bins ...
TH2F * fhConvDeltaEtaPhiMCString
! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is string ...
const Double_t ptmin
TH2F * fhConvDistEn
! Approx distance to vertex vs Energy
Float_t fMassCutTight
Mass cut for the conversion pairs selection, tighter.
virtual Double_t GetEventWeight() const
Bool_t fFillClusterConvDistHisto
Fill histograms with calculated conversion distance with data clusters.
TH2F * fhConvDistEtaCutAll
! Approx distance to vertex vs rapidity, dEta < fConvDEtaCut, M < 20 MeV/c^2, A < fConvAsymCut ...
virtual TObjArray * GetPHOSClusters() const
Float_t GetHistoEtaMin() const
TH2F * fhConvDistMCConversion
! Calculated conversion distance vs real distance to vertex
virtual void SetLabel(Int_t l)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
virtual void SetTag(Int_t tag)
TH2F * fhConvDeltaEtaMCAntiProton
! Small mass cluster pairs, correlation in eta, origin of both clusters is anti proton ...
TH2F * fhConvDistPhiCutAll
! Approx distance to vertex vs azimuth, dEta < fConvDEtaCut, M < 20 MeV/c^2, A < fConvAsymCut ...
TH1F * fhPtAntiProtonTagged
! Number of identified gamma from AntiProtons gamma, tagged as conversion
TH2F * fhConvPtMCAntiProton
! Small mass cluster pairs, pt of pairs, origin of both clusters is anti proton
TH2F * fhConvDistEtaCutMass
! Approx distance to vertex vs rapidity, M < 20 MeV/c^2
TH2F * fhConvAsymMCAntiProton
! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is anti proton ...
virtual void SetCaloLabel(Int_t a, Int_t b)
Float_t GetHistoEtaMax() const
Int_t GetHistoPtBins() const
virtual void SetTagged(Bool_t tag)
TH2F * fhConvM02MCString
! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is string ...
Bool_t fAddConvertedPairsToAOD
Put Converted pairs in AOD.
TLorentzVector fMomentum
! Cluster momentum
TH2F * fhConvDistMassCutAsy
! Approx distance to vertex vs mass, A < fConvAsymCut
TH2F * fhConvAsymMCAntiNeutron
! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is anti neutron ...
TH2F * fhConvDistPhiCutEta
! Approx distance to vertex vs azimuth
virtual Int_t GetModuleNumber(AliCaloTrackParticle *part) const
const Double_t etamax
virtual Double_t Phi() const
Conversions pairs clusters analysis.
virtual AliMCAnalysisUtils * GetMCAnalysisUtils()
TVector3 fProdVertex
! Production vertex
Float_t fConvAsymCut
Select conversion pairs when asymmetry is smaller than cut.
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
TH2F * fhConvAsymMCConversion
! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is conversion ...
const char Option_t
Definition: External.C:48
Float_t fMassCut
Mass cut for the conversion pairs selection.
Int_t CheckCommonAncestor(Int_t index1, Int_t index2, const AliMCEvent *mcevent, Int_t &ancPDG, Int_t &ancStatus, TLorentzVector &momentum, TVector3 &prodVertex)
virtual AliVCluster * FindCluster(TObjArray *clusters, Int_t clId, Int_t &iclus, Int_t first=0)
virtual Float_t GetM02() const
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Float_t GetHistoPhiMax() const
virtual AliCaloTrackReader * GetReader() const
bool Bool_t
Definition: External.C:53
Int_t GetHistoEtaBins() const
TH2F * fhConvDistEnCutMass
! Approx distance to vertex vs Energy, M < 20 MeV/c^2
virtual Int_t GetTag() const
virtual TObjArray * GetEMCALClusters() const
virtual Int_t GetIdentifiedParticleType() const
Int_t nptbins
virtual Int_t GetCaloLabel(Int_t i) const
TH2F * fhConvDeltaEtaPhiMCAntiNeutron
! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is anti neutron ...
TH2F * fhConvDeltaPhiMCString
! Small mass cluster pairs, correlation in phi, origin of both clusters is string ...
const Double_t phimin
TH2F * fhConvDistPhiCutAsy
! Approx distance to vertex vs azimuth, A < fConvAsymCut