1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
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
24 //-- Author: Gustavo Conesa (LPSC-IN2P3-CNRS)
25 //////////////////////////////////////////////////////////////////////////////
28 // --- ROOT system ---
31 #include <TClonesArray.h>
32 #include <TObjString.h>
33 #include "TParticle.h"
34 #include "TDatabasePDG.h"
36 // --- Analysis system ---
37 #include "AliAnaPhotonConvInCalo.h"
38 #include "AliCaloTrackReader.h"
40 #include "AliCaloPID.h"
41 #include "AliMCAnalysisUtils.h"
42 #include "AliFiducialCut.h"
43 #include "AliVCluster.h"
44 #include "AliAODMCParticle.h"
46 ClassImp(AliAnaPhotonConvInCalo)
48 //________________________________________________
49 AliAnaPhotonConvInCalo::AliAnaPhotonConvInCalo() :
50 AliAnaCaloTrackCorrBaseClass(),
51 fRemoveConvertedPair(kFALSE),
52 fAddConvertedPairsToAOD(kFALSE),
54 fConvAsymCut(1.), fConvDEtaCut(2.),
55 fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
56 fMomentum(), fProdVertex(),
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),
67 fhPtConversionTagged(0), fhPtAntiNeutronTagged(0),
68 fhPtAntiProtonTagged(0), fhPtUnknownTagged(0),
70 fhConvDeltaEtaMCConversion(0), fhConvDeltaPhiMCConversion(0), fhConvDeltaEtaPhiMCConversion(0),
71 fhConvAsymMCConversion(0), fhConvPtMCConversion(0),
72 fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
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)
87 //Initialize parameters
92 //_____________________________________________________
93 TObjString * AliAnaPhotonConvInCalo::GetAnalysisCuts()
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] ;
100 snprintf(onePar,buffersize,"--- AliAnaPhotonConvInCalo---:") ;
102 snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)",
103 fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
106 return new TObjString(parList) ;
109 //_______________________________________________________
110 TList * AliAnaPhotonConvInCalo::GetCreateOutputObjects()
112 // Create histograms to be saved in output file and
113 // store them in outputContainer
114 TList * outputContainer = new TList() ;
115 outputContainer->SetName("PhotonConvInCaloHistos") ;
117 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
118 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
119 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
434 return outputContainer ;
438 //___________________________________________
439 void AliAnaPhotonConvInCalo::InitParameters()
441 //Initialize the parameters of the analysis.
442 AddToHistogramsName("AnaPhotonConvInCalo_");
444 fMassCut = 0.03; //30 MeV
445 fRemoveConvertedPair = kFALSE;
446 fAddConvertedPairsToAOD = kFALSE;
450 //_________________________________________________
451 void AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()
453 //Do conversion photon analysis and fill aods
455 //Loop on stored AOD photons
456 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
457 AliDebug(1,Form("AOD branch entries %d", naod));
459 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
460 Bool_t * indexConverted = new Bool_t[naod];
461 for (Int_t i = 0; i < naod; i++) indexConverted[i] = kFALSE;
463 for(Int_t iaod = 0; iaod < naod ; iaod++)
465 AliAODPWG4Particle* calo = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
467 Bool_t bConverted = kFALSE;
470 //Check if set previously as converted couple, if so skip its use.
471 if (indexConverted[iaod]) continue;
473 // Second cluster loop
474 AliAODPWG4Particle* calo2 = 0;
475 for(Int_t jaod = iaod + 1 ; jaod < naod ; jaod++)
477 //Check if set previously as converted couple, if so skip its use.
478 if (indexConverted[jaod]) continue;
479 //printf("Check Conversion indeces %d and %d\n",iaod,jaod);
480 calo2 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(jaod));
482 //................................................
483 //Get mass of pair, if small, take this pair.
484 Float_t pairM = calo->GetPairMass(calo2);
485 //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
488 calo->SetTagged(kFALSE);
489 id2 = calo2->GetCaloLabel(0);
490 Float_t asymmetry = TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E());
491 Float_t dPhi = (calo->Momentum())->Phi()-(calo2->Momentum())->Phi();
492 Float_t dEta = (calo->Momentum())->Eta()-(calo2->Momentum())->Eta();
494 //...............................................
495 //Fill few histograms with kinematics of the pair
496 //FIXME, move all this to MakeAnalysisFillHistograms ...
498 fhConvDeltaEta ->Fill( pairM, dPhi );
499 fhConvDeltaPhi ->Fill( pairM, dEta );
500 fhConvAsym ->Fill( pairM, asymmetry );
501 fhConvDeltaEtaPhi->Fill( dEta , dPhi );
502 fhConvPt ->Fill( pairM, (calo->Momentum())->Pt()+(calo2->Momentum())->Pt());
504 //Estimate conversion distance, T. Awes, M. Ivanov
505 //Under the assumption that the pair has zero mass, and that each electron
506 //of the pair has the same momentum, they will each have the same bend radius
507 //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
508 //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
509 //R = E/1.5 [cm]. Under these assumptions, the distance from the conversion
510 //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
511 //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
512 //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
515 TObjArray * clusters = 0;
516 if(calo->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
517 else clusters = GetPHOSClusters();
521 AliVCluster *cluster1 = FindCluster(clusters,calo ->GetCaloLabel(0),iclus);
522 AliVCluster *cluster2 = FindCluster(clusters,calo2->GetCaloLabel(0),iclus);
525 cluster1->GetPosition(pos1);
527 cluster2->GetPosition(pos2);
528 Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
529 (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
530 (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
532 Float_t convDist = TMath::Sqrt(calo->E() *clustDist*0.01/0.15);
533 Float_t convDist2 = TMath::Sqrt(calo2->E()*clustDist*0.01/0.15);
534 //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,calo->E(),convDist,calo2->E(),convDist2);
535 AliDebug(2,Form("Pair with mass %2.3f < %2.3f, %1.2f < dPhi %2.2f < %2.2f, dEta %f < %2.2f, asymmetry %2.2f< %2.2f; \n cluster1 id %d, e %2.3f SM %d, eta %2.3f, phi %2.3f ; \n cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f",
536 pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
537 calo->GetCaloLabel(0),calo->E(),GetCaloUtils()->GetModuleNumber(calo,GetReader()->GetInputEvent()), calo->Eta(), calo->Phi(),
538 id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2,GetReader()->GetInputEvent()),calo2->Eta(), calo2->Phi()));
540 fhConvDistEta ->Fill(calo ->Eta(),convDist );
541 fhConvDistEta ->Fill(calo2->Eta(),convDist2);
542 fhConvDistEn ->Fill(calo ->E(), convDist );
543 fhConvDistEn ->Fill(calo2->E(), convDist2);
544 fhConvDistMass->Fill(pairM, convDist );
548 fhConvDistEtaCutEta ->Fill(calo->Eta(), convDist );
549 fhConvDistEtaCutEta ->Fill(calo2->Eta(),convDist2);
550 fhConvDistEnCutEta ->Fill(calo->E(), convDist );
551 fhConvDistEnCutEta ->Fill(calo2->E(), convDist2);
552 fhConvDistMassCutEta->Fill(pairM, convDist );
556 fhConvDistEtaCutMass ->Fill(calo ->Eta(), convDist );
557 fhConvDistEtaCutMass ->Fill(calo2->Eta(), convDist2);
558 fhConvDistEnCutMass ->Fill(calo ->E(), convDist );
559 fhConvDistEnCutMass ->Fill(calo2->E(), convDist2);
563 fhConvDistEtaCutAsy ->Fill(calo ->Eta(), convDist );
564 fhConvDistEtaCutAsy ->Fill(calo2->Eta(), convDist2);
565 fhConvDistEnCutAsy ->Fill(calo ->E(), convDist );
566 fhConvDistEnCutAsy ->Fill(calo2->E(), convDist2);
571 //...............................................
572 //Select pairs in a eta-phi window
573 if(TMath::Abs(dEta) < fConvDEtaCut &&
574 TMath::Abs(dPhi) < fConvDPhiMaxCut &&
575 TMath::Abs(dPhi) > fConvDPhiMinCut &&
576 asymmetry < fConvAsymCut )
578 indexConverted[iaod] = kTRUE;
579 indexConverted[jaod] = kTRUE;
582 //printf("Accepted? %d\n",bConverted);
583 //...........................................
584 //Fill more histograms, simulated data
585 //FIXME, move all this to MakeAnalysisFillHistograms ...
588 //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
591 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster1->GetLabel(), cluster2->GetLabel(),
592 GetReader(), ancPDG, ancStatus, fMomentum, fProdVertex);
594 // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
595 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
597 Int_t tag1 = calo ->GetTag();
598 Int_t tag2 = calo2->GetTag();
599 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion))
601 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1)
603 fhConvDeltaEtaMCConversion ->Fill( pairM, dEta );
604 fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi );
605 fhConvAsymMCConversion ->Fill( pairM, asymmetry );
606 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi );
607 fhConvPtMCConversion ->Fill( pairM, calo->Pt()+calo2->Pt());
608 fhConvDispersionMCConversion ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
609 fhConvM02MCConversion ->Fill( cluster1->GetM02(), cluster2->GetM02());
610 fhConvDistMCConversion ->Fill( convDist , fProdVertex.Mag() );
611 fhConvDistMCConversion ->Fill( convDist2, fProdVertex.Mag() );
613 if(dEta<0.05 && pairM<0.01 && asymmetry<0.1)
615 fhConvDistMCConversionCuts->Fill( convDist , fProdVertex.Mag() );
616 fhConvDistMCConversionCuts->Fill( convDist2, fProdVertex.Mag() );
621 else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiNeutron))
623 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1)
625 fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta );
626 fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi );
627 fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry );
628 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi );
629 fhConvPtMCAntiNeutron ->Fill( pairM, calo->Pt()+calo2->Pt());
630 fhConvDispersionMCAntiNeutron ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
631 fhConvM02MCAntiNeutron ->Fill( cluster1->GetM02(), cluster2->GetM02());
634 else if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCAntiProton))
636 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1)
638 fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta );
639 fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi );
640 fhConvAsymMCAntiProton ->Fill( pairM, asymmetry );
641 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi );
642 fhConvPtMCAntiProton ->Fill( pairM, calo->Pt()+calo2->Pt());
643 fhConvDispersionMCAntiProton ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
644 fhConvM02MCAntiProton ->Fill( cluster1->GetM02(), cluster2->GetM02());
648 //Pairs coming from fragmenting pairs.
649 if( ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) )
651 fhConvDeltaEtaMCString ->Fill( pairM, dPhi);
652 fhConvDeltaPhiMCString ->Fill( pairM, dPhi);
653 fhConvAsymMCString ->Fill( pairM, TMath::Abs(calo->E()-calo2->E())/(calo->E()+calo2->E()) );
654 fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi);
655 fhConvPtMCString ->Fill( pairM, calo->Pt()+calo2->Pt());
656 fhConvDispersionMCString ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
657 fhConvM02MCString ->Fill( cluster1->GetM02(), cluster2->GetM02());
667 //..........................................................................................................
668 //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
672 if(fAddConvertedPairsToAOD)
674 //Create AOD of pair analysis
675 fMomentum = *(calo->Momentum())+*(calo2->Momentum());
676 AliAODPWG4Particle aodpair = AliAODPWG4Particle(fMomentum);
677 aodpair.SetLabel(calo->GetLabel());
679 //printf("Index %d, Id %d\n",iaod, calo->GetID());
680 //Set the indeces of the original caloclusters
681 aodpair.SetCaloLabel(calo->GetCaloLabel(0),id2);
682 aodpair.SetDetectorTag(calo->GetDetectorTag());
683 aodpair.SetIdentifiedParticleType(calo->GetIdentifiedParticleType());
684 aodpair.SetTag(calo ->GetTag());
685 aodpair.SetTagged(kTRUE);
686 //Add AOD with pair object to aod branch
687 AddAODParticle(aodpair);
688 //printf("\t \t both added pair\n");
691 //Do not add the current calocluster
692 if(!fRemoveConvertedPair)
694 //printf("TAGGED\n");
695 //Tag this cluster as likely conversion
696 calo->SetTagged(kTRUE);
702 // Remove entries identified as conversion electrons
703 // Revise if this is OK
704 if(fRemoveConvertedPair || fAddConvertedPairsToAOD)
706 for(Int_t iaod = 0; iaod < naod ; iaod++)
708 if(indexConverted[iaod])GetOutputAODBranch()->RemoveAt(iaod);
710 GetOutputAODBranch()->Compress();
713 delete [] indexConverted;
715 AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
719 //________________________________________________________
720 void AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms()
724 //-------------------------------------------------------------------
725 // Access MC information in stack if requested, check that it exists.
726 AliStack * stack = 0x0;
727 TParticle * primary = 0x0;
728 TClonesArray * mcparticles = 0x0;
729 AliAODMCParticle * aodprimary = 0x0;
733 if(GetReader()->ReadStack())
735 stack = GetMCStack() ;
738 AliFatal("Stack not available, is the MC handler called? STOP");
742 else if(GetReader()->ReadAODMCParticles())
744 //Get the list of MC particles
745 mcparticles = GetReader()->GetAODMCParticles();
748 AliFatal("Standard MCParticles not available!");
754 //----------------------------------
755 //Loop on stored AOD photons
756 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
757 AliDebug(1,Form("AOD branch entries %d", naod));
759 for(Int_t iaod = 0; iaod < naod ; iaod++)
761 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
762 //Int_t pdg = ph->GetIdentifiedParticleType();
766 AliDebug(2,Form("ID Photon: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta())) ;
767 //................................
768 //Fill photon histograms
769 Float_t ptcluster = ph->Pt();
770 Float_t phicluster = ph->Phi();
771 Float_t etacluster = ph->Eta();
772 Float_t ecluster = ph->E();
774 fhPtPhotonConv->Fill(ptcluster);
775 if(ecluster > 0.5) fhEtaPhiPhotonConv ->Fill(etacluster, phicluster);
776 else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
778 //.......................................
779 //Play with the MC data if available
782 //....................................................................
783 // Access MC information in stack if requested, check that it exists.
784 Int_t label =ph->GetLabel();
787 AliDebug(1,Form("*** bad label ***: label %d", label));
792 //Float_t ptprim = 0;
793 if(GetReader()->ReadStack())
795 if(label >= stack->GetNtrack())
797 AliDebug(1,Form("*** large label ***: label %d, n tracks %d", label, stack->GetNtrack()));
801 primary = stack->Particle(label);
804 AliDebug(1,Form("*** no primary ***: label %d", label));
807 //eprim = primary->Energy();
808 //ptprim = primary->Pt();
811 else if(GetReader()->ReadAODMCParticles())
813 if(label >= mcparticles->GetEntriesFast())
815 AliDebug(2,Form("*** large label ***: label %d, n tracks %d",label, mcparticles->GetEntriesFast()));
820 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
824 AliDebug(2,Form("*** no primary ***: label %d", label));
828 //eprim = aodprimary->E();
829 //ptprim = aodprimary->Pt();
833 Int_t tag =ph->GetTag();
835 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
837 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
839 fhPtConversionTagged ->Fill(ptcluster);
842 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
844 fhPtAntiNeutronTagged ->Fill(ptcluster);
846 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
848 fhPtAntiProtonTagged ->Fill(ptcluster);
852 fhPtUnknownTagged ->Fill(ptcluster);
855 }//Histograms with MC
856 }// tagged by conversion
862 //________________________________________________________
863 void AliAnaPhotonConvInCalo::Print(const Option_t * opt) const
865 //Print some relevant parameters set for the analysis
870 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
871 AliAnaCaloTrackCorrBaseClass::Print(" ");
873 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
874 printf("Conversion pair mass cut = %f\n",fMassCut);
875 printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
876 fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);