]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPhoton.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Reader: Add option to remove or not event with primary vertex not reconstructed
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPhoton.cxx
CommitLineData
a3aebfff 1 /**************************************************************************
1c5acb87 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 *
cadbb0f3 9 * without fee, provided that the above copyright notice appears in all *
1c5acb87 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/* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18//
19// Class for the photon identification.
20// Clusters from calorimeters are identified as photons
21// and kept in the AOD. Few histograms produced.
6175da48 22// Produces input for other analysis classes like AliAnaPi0,
23// AliAnaParticleHadronCorrelation ...
1c5acb87 24//
25// -- Author: Gustavo Conesa (LNF-INFN)
26//////////////////////////////////////////////////////////////////////////////
27
28
29// --- ROOT system ---
30#include <TH2F.h>
2244659d 31#include <TH3D.h>
477d6cee 32#include <TClonesArray.h>
0c1383b5 33#include <TObjString.h>
477d6cee 34//#include <Riostream.h>
123fc3bd 35#include "TParticle.h"
6175da48 36#include "TDatabasePDG.h"
1c5acb87 37
38// --- Analysis system ---
39#include "AliAnaPhoton.h"
40#include "AliCaloTrackReader.h"
123fc3bd 41#include "AliStack.h"
1c5acb87 42#include "AliCaloPID.h"
6639984f 43#include "AliMCAnalysisUtils.h"
ff45398a 44#include "AliFiducialCut.h"
0ae57829 45#include "AliVCluster.h"
591cc579 46#include "AliAODMCParticle.h"
c8fe2783 47#include "AliMixedEvent.h"
48
1c5acb87 49
50ClassImp(AliAnaPhoton)
51
52//____________________________________________________________________________
53 AliAnaPhoton::AliAnaPhoton() :
54 AliAnaPartCorrBaseClass(), fCalorimeter(""),
a3aebfff 55 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
2ac125bf 56 fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
20218aea 57 fCheckConversion(kFALSE), fRemoveConvertedPair(kFALSE), fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
41121cfe 58 fConvAsymCut(1.), fConvDEtaCut(2.),fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
6175da48 59 //fhVertex(0),
60 fhNtraNclu(0), fhNCellsPt(0),
20218aea 61 fhEPhoton(0), fhPtPhoton(0), fhPhiPhoton(0), fhEtaPhoton(0), fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
6175da48 62 fhPtPhotonConv(0), fhEtaPhiPhotonConv(0),fhEtaPhi05PhotonConv(0),
41121cfe 63 fhConvDeltaEta(0), fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0), fhConvAsym(0), fhConvPt(0),
7175a03a 64 //MC
123fc3bd 65 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
0ae57829 66 fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
1c5acb87 67 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
68 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
6639984f 69 fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
1c5acb87 70 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
6175da48 71 fhPtOtherDecay(0), fhPhiOtherDecay(0), fhEtaOtherDecay(0),
41121cfe 72 fhPtConversion(0), fhPhiConversion(0), fhEtaConversion(0),fhEtaPhiConversion(0),fhEtaPhi05Conversion(0),
6175da48 73 fhPtAntiNeutron(0), fhPhiAntiNeutron(0), fhEtaAntiNeutron(0),
74 fhPtAntiProton(0), fhPhiAntiProton(0), fhEtaAntiProton(0),
75 fhPtUnknown(0), fhPhiUnknown(0), fhEtaUnknown(0),
41121cfe 76 fhPtConversionTagged(0), fhPtAntiNeutronTagged(0), fhPtAntiProtonTagged(0), fhPtUnknownTagged(0),
6175da48 77 fhConvDeltaEtaMCConversion(0), fhConvDeltaPhiMCConversion(0), fhConvDeltaEtaPhiMCConversion(0), fhConvAsymMCConversion(0), fhConvPtMCConversion(0), fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
78 fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0), fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0), fhConvDispersionMCAntiNeutron(0),fhConvM02MCAntiNeutron(0),
79 fhConvDeltaEtaMCAntiProton(0), fhConvDeltaPhiMCAntiProton(0), fhConvDeltaEtaPhiMCAntiProton(0), fhConvAsymMCAntiProton(0), fhConvPtMCAntiProton(0), fhConvDispersionMCAntiProton(0), fhConvM02MCAntiProton(0),
80 fhConvDeltaEtaMCString(0), fhConvDeltaPhiMCString(0), fhConvDeltaEtaPhiMCString(0), fhConvAsymMCString(0), fhConvPtMCString(0), fhConvDispersionMCString(0), fhConvM02MCString(0)
1c5acb87 81{
82 //default ctor
83
84 //Initialize parameters
85 InitParameters();
86
5ae09196 87}//____________________________________________________________________________
1c5acb87 88AliAnaPhoton::~AliAnaPhoton()
89{
90 //dtor
91
92}
93
0c1383b5 94//________________________________________________________________________
95TObjString * AliAnaPhoton::GetAnalysisCuts()
96{
97 //Save parameters used for analysis
98 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 99 const Int_t buffersize = 255;
100 char onePar[buffersize] ;
0c1383b5 101
5ae09196 102 snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
0c1383b5 103 parList+=onePar ;
5ae09196 104 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
0c1383b5 105 parList+=onePar ;
5ae09196 106 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
0c1383b5 107 parList+=onePar ;
5ae09196 108 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
0c1383b5 109 parList+=onePar ;
5ae09196 110 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
0c1383b5 111 parList+=onePar ;
5ae09196 112 snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
0c1383b5 113 parList+=onePar ;
41121cfe 114 snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
115 fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
6175da48 116 parList+=onePar ;
0c1383b5 117
118 //Get parameters set in base class.
119 parList += GetBaseParametersList() ;
120
121 //Get parameters set in PID class.
122 parList += GetCaloPID()->GetPIDParametersList() ;
123
124 //Get parameters set in FiducialCut class (not available yet)
125 //parlist += GetFidCut()->GetFidCutParametersList()
126
127 return new TObjString(parList) ;
128}
129
1c5acb87 130
131//________________________________________________________________________
132TList * AliAnaPhoton::GetCreateOutputObjects()
133{
477d6cee 134 // Create histograms to be saved in output file and
135 // store them in outputContainer
136 TList * outputContainer = new TList() ;
137 outputContainer->SetName("PhotonHistos") ;
4a745797 138
5a2dbc3c 139 Int_t nptbins = GetHistoPtBins();
140 Int_t nphibins = GetHistoPhiBins();
141 Int_t netabins = GetHistoEtaBins();
477d6cee 142 Float_t ptmax = GetHistoPtMax();
143 Float_t phimax = GetHistoPhiMax();
144 Float_t etamax = GetHistoEtaMax();
145 Float_t ptmin = GetHistoPtMin();
146 Float_t phimin = GetHistoPhiMin();
147 Float_t etamin = GetHistoEtaMin();
148
149 //Histograms of highest Photon identified in Event
6175da48 150// fhVertex = new TH3D ("Vertex","vertex position", 20,-10.,10., 20,-10.,10., 80,-40.,40.);
151// fhVertex->SetXTitle("X");
152// fhVertex->SetYTitle("Y");
153// fhVertex->SetZTitle("Z");
154// outputContainer->Add(fhVertex);
2244659d 155
156 fhNtraNclu = new TH2F ("hNtracksNcluster","# of tracks vs # of clusters", 500,0,500, 500,0,500);
157 fhNtraNclu->SetXTitle("# of tracks");
158 fhNtraNclu->SetYTitle("# of clusters");
159 outputContainer->Add(fhNtraNclu);
160
ba1eeb1f 161 fhNCellsPt = new TH2F ("hNCellsPt","# of cells in cluster vs E of clusters", nptbins,ptmin, ptmax, 10,0,10);
162 fhNCellsPt->SetXTitle("p_{T} (GeV/c)");
6175da48 163 fhNCellsPt->SetYTitle("# of cells in cluster");
164 outputContainer->Add(fhNCellsPt);
165
20218aea 166 fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
167 fhEPhoton->SetYTitle("N");
168 fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
169 outputContainer->Add(fhEPhoton) ;
170
171 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
477d6cee 172 fhPtPhoton->SetYTitle("N");
173 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
174 outputContainer->Add(fhPtPhoton) ;
175
176 fhPhiPhoton = new TH2F
20218aea 177 ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
6175da48 178 fhPhiPhoton->SetYTitle("#phi (rad)");
477d6cee 179 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
180 outputContainer->Add(fhPhiPhoton) ;
181
182 fhEtaPhoton = new TH2F
20218aea 183 ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 184 fhEtaPhoton->SetYTitle("#eta");
185 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
186 outputContainer->Add(fhEtaPhoton) ;
187
6175da48 188 fhEtaPhiPhoton = new TH2F
189 ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
190 fhEtaPhiPhoton->SetYTitle("#phi (rad)");
191 fhEtaPhiPhoton->SetXTitle("#eta");
192 outputContainer->Add(fhEtaPhiPhoton) ;
20218aea 193 if(GetMinPt() < 0.5){
194 fhEtaPhi05Photon = new TH2F
195 ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
196 fhEtaPhi05Photon->SetYTitle("#phi (rad)");
197 fhEtaPhi05Photon->SetXTitle("#eta");
198 outputContainer->Add(fhEtaPhi05Photon) ;
199 }
6175da48 200
201 //Conversion
20218aea 202 if(fCheckConversion){
203 fhPtPhotonConv = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
204 fhPtPhotonConv->SetYTitle("N");
205 fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
206 outputContainer->Add(fhPtPhotonConv) ;
207
208 fhEtaPhiPhotonConv = new TH2F
209 ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
210 fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
211 fhEtaPhiPhotonConv->SetXTitle("#eta");
212 outputContainer->Add(fhEtaPhiPhotonConv) ;
213 if(GetMinPt() < 0.5){
214 fhEtaPhi05PhotonConv = new TH2F
215 ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
216 fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
217 fhEtaPhi05PhotonConv->SetXTitle("#eta");
218 outputContainer->Add(fhEtaPhi05PhotonConv) ;
219 }
220
221 fhConvDeltaEta = new TH2F
222 ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
223 fhConvDeltaEta->SetYTitle("#Delta #eta");
224 fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
225 outputContainer->Add(fhConvDeltaEta) ;
226
227 fhConvDeltaPhi = new TH2F
228 ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
229 fhConvDeltaPhi->SetYTitle("#Delta #phi");
230 fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
231 outputContainer->Add(fhConvDeltaPhi) ;
232
233 fhConvDeltaEtaPhi = new TH2F
234 ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
235 fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
236 fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
237 outputContainer->Add(fhConvDeltaEtaPhi) ;
238
239 fhConvAsym = new TH2F
240 ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
241 fhConvAsym->SetYTitle("Asymmetry");
242 fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
243 outputContainer->Add(fhConvAsym) ;
244
245 fhConvPt = new TH2F
246 ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
247 fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
248 fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
249 outputContainer->Add(fhConvPt) ;
250 }
6175da48 251
477d6cee 252 if(IsDataMC()){
123fc3bd 253 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
254 fhDeltaE->SetXTitle("#Delta E (GeV)");
255 outputContainer->Add(fhDeltaE);
256
257 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
258 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
259 outputContainer->Add(fhDeltaPt);
260
261 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
262 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
263 outputContainer->Add(fhRatioE);
477d6cee 264
123fc3bd 265 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
266 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
267 outputContainer->Add(fhRatioPt);
268
269 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 270 fh2E->SetXTitle("E_{rec} (GeV)");
271 fh2E->SetYTitle("E_{gen} (GeV)");
123fc3bd 272 outputContainer->Add(fh2E);
273
274 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 275 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
276 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
123fc3bd 277 outputContainer->Add(fh2Pt);
278
c8fe2783 279 fhPtMCPhoton = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
280 fhPtMCPhoton->SetYTitle("N");
281 fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
282 outputContainer->Add(fhPtMCPhoton) ;
283
284 fhPhiMCPhoton = new TH2F
5ae09196 285 ("hPhiMCPhoton","#phi_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
c8fe2783 286 fhPhiMCPhoton->SetYTitle("#phi");
287 fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
288 outputContainer->Add(fhPhiMCPhoton) ;
289
290 fhEtaMCPhoton = new TH2F
5ae09196 291 ("hEtaMCPhoton","#eta_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
c8fe2783 292 fhEtaMCPhoton->SetYTitle("#eta");
293 fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
294 outputContainer->Add(fhEtaMCPhoton) ;
295
591cc579 296 fhPtPrompt = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 297 fhPtPrompt->SetYTitle("N");
298 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
299 outputContainer->Add(fhPtPrompt) ;
300
301 fhPhiPrompt = new TH2F
5ae09196 302 ("hPhiMCPrompt","#phi_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 303 fhPhiPrompt->SetYTitle("#phi");
304 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
305 outputContainer->Add(fhPhiPrompt) ;
306
307 fhEtaPrompt = new TH2F
5ae09196 308 ("hEtaMCPrompt","#eta_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 309 fhEtaPrompt->SetYTitle("#eta");
310 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
311 outputContainer->Add(fhEtaPrompt) ;
312
591cc579 313 fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 314 fhPtFragmentation->SetYTitle("N");
315 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
316 outputContainer->Add(fhPtFragmentation) ;
317
318 fhPhiFragmentation = new TH2F
5ae09196 319 ("hPhiMCFragmentation","#phi_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 320 fhPhiFragmentation->SetYTitle("#phi");
321 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
322 outputContainer->Add(fhPhiFragmentation) ;
323
324 fhEtaFragmentation = new TH2F
5ae09196 325 ("hEtaMCFragmentation","#eta_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 326 fhEtaFragmentation->SetYTitle("#eta");
327 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
328 outputContainer->Add(fhEtaFragmentation) ;
329
a3aebfff 330 fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 331 fhPtISR->SetYTitle("N");
332 fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
333 outputContainer->Add(fhPtISR) ;
334
335 fhPhiISR = new TH2F
a3aebfff 336 ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 337 fhPhiISR->SetYTitle("#phi");
338 fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
339 outputContainer->Add(fhPhiISR) ;
340
341 fhEtaISR = new TH2F
5ae09196 342 ("hEtaMCISR","#eta_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 343 fhEtaISR->SetYTitle("#eta");
344 fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
345 outputContainer->Add(fhEtaISR) ;
346
591cc579 347 fhPtPi0Decay = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 348 fhPtPi0Decay->SetYTitle("N");
349 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
350 outputContainer->Add(fhPtPi0Decay) ;
351
352 fhPhiPi0Decay = new TH2F
5ae09196 353 ("hPhiMCPi0Decay","#phi_{#gamma}, #pi^{0} decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 354 fhPhiPi0Decay->SetYTitle("#phi");
355 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
356 outputContainer->Add(fhPhiPi0Decay) ;
357
358 fhEtaPi0Decay = new TH2F
5ae09196 359 ("hEtaMCPi0Decay","#eta_{#gamma}, #pi^{0} #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 360 fhEtaPi0Decay->SetYTitle("#eta");
361 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
362 outputContainer->Add(fhEtaPi0Decay) ;
363
a3aebfff 364 fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 365 fhPtOtherDecay->SetYTitle("N");
366 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
367 outputContainer->Add(fhPtOtherDecay) ;
368
369 fhPhiOtherDecay = new TH2F
5ae09196 370 ("hPhiMCOtherDecay","#phi_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 371 fhPhiOtherDecay->SetYTitle("#phi");
372 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
373 outputContainer->Add(fhPhiOtherDecay) ;
374
375 fhEtaOtherDecay = new TH2F
5ae09196 376 ("hEtaMCOtherDecay","#eta_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 377 fhEtaOtherDecay->SetYTitle("#eta");
378 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
379 outputContainer->Add(fhEtaOtherDecay) ;
380
a3aebfff 381 fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 382 fhPtConversion->SetYTitle("N");
383 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
384 outputContainer->Add(fhPtConversion) ;
385
386 fhPhiConversion = new TH2F
5ae09196 387 ("hPhiMCConversion","#phi_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 388 fhPhiConversion->SetYTitle("#phi");
389 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
390 outputContainer->Add(fhPhiConversion) ;
391
392 fhEtaConversion = new TH2F
5ae09196 393 ("hEtaMCConversion","#eta_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 394 fhEtaConversion->SetYTitle("#eta");
395 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
396 outputContainer->Add(fhEtaConversion) ;
397
6175da48 398 fhEtaPhiConversion = new TH2F
399 ("hEtaPhiConversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
400 fhEtaPhiConversion->SetYTitle("#phi (rad)");
401 fhEtaPhiConversion->SetXTitle("#eta");
402 outputContainer->Add(fhEtaPhiConversion) ;
403
404 fhEtaPhi05Conversion = new TH2F
405 ("hEtaPhi05Conversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
406 fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
407 fhEtaPhi05Conversion->SetXTitle("#eta");
408 outputContainer->Add(fhEtaPhi05Conversion) ;
409
410 fhPtAntiNeutron = new TH1F("hPtMCAntiNeutron","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
411 fhPtAntiNeutron->SetYTitle("N");
412 fhPtAntiNeutron->SetXTitle("p_{T #gamma}(GeV/c)");
413 outputContainer->Add(fhPtAntiNeutron) ;
414
415 fhPhiAntiNeutron = new TH2F
416 ("hPhiMCAntiNeutron","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
417 fhPhiAntiNeutron->SetYTitle("#phi");
418 fhPhiAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
419 outputContainer->Add(fhPhiAntiNeutron) ;
420
421 fhEtaAntiNeutron = new TH2F
422 ("hEtaMCAntiNeutron","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
423 fhEtaAntiNeutron->SetYTitle("#eta");
424 fhEtaAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
425 outputContainer->Add(fhEtaAntiNeutron) ;
426
41121cfe 427 fhPtAntiProton = new TH1F("hPtMCAntiProton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
6175da48 428 fhPtAntiProton->SetYTitle("N");
429 fhPtAntiProton->SetXTitle("p_{T #gamma}(GeV/c)");
430 outputContainer->Add(fhPtAntiProton) ;
431
432 fhPhiAntiProton = new TH2F
433 ("hPhiMCAntiProton","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
434 fhPhiAntiProton->SetYTitle("#phi");
435 fhPhiAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
436 outputContainer->Add(fhPhiAntiProton) ;
437
438 fhEtaAntiProton = new TH2F
439 ("hEtaMCAntiProton","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
440 fhEtaAntiProton->SetYTitle("#eta");
441 fhEtaAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
442 outputContainer->Add(fhEtaAntiProton) ;
443
a3aebfff 444 fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 445 fhPtUnknown->SetYTitle("N");
446 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
447 outputContainer->Add(fhPtUnknown) ;
448
449 fhPhiUnknown = new TH2F
5ae09196 450 ("hPhiMCUnknown","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 451 fhPhiUnknown->SetYTitle("#phi");
452 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
453 outputContainer->Add(fhPhiUnknown) ;
454
455 fhEtaUnknown = new TH2F
5ae09196 456 ("hEtaMCUnknown","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 457 fhEtaUnknown->SetYTitle("#eta");
458 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
459 outputContainer->Add(fhEtaUnknown) ;
a3aebfff 460
20218aea 461 if(fCheckConversion){
462 fhPtConversionTagged = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
463 fhPtConversionTagged->SetYTitle("N");
464 fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
465 outputContainer->Add(fhPtConversionTagged) ;
466
467
468 fhPtAntiNeutronTagged = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
469 fhPtAntiNeutronTagged->SetYTitle("N");
470 fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
471 outputContainer->Add(fhPtAntiNeutronTagged) ;
472
473 fhPtAntiProtonTagged = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
474 fhPtAntiProtonTagged->SetYTitle("N");
475 fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
476 outputContainer->Add(fhPtAntiProtonTagged) ;
477
478 fhPtUnknownTagged = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
479 fhPtUnknownTagged->SetYTitle("N");
480 fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
481 outputContainer->Add(fhPtUnknownTagged) ;
482
483 fhConvDeltaEtaMCConversion = new TH2F
484 ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
485 fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
486 fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
487 outputContainer->Add(fhConvDeltaEtaMCConversion) ;
488
489 fhConvDeltaPhiMCConversion = new TH2F
490 ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
491 fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
492 fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
493 outputContainer->Add(fhConvDeltaPhiMCConversion) ;
494
495 fhConvDeltaEtaPhiMCConversion = new TH2F
496 ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
497 fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
498 fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
499 outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
500
501 fhConvAsymMCConversion = new TH2F
502 ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
503 fhConvAsymMCConversion->SetYTitle("Asymmetry");
504 fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
505 outputContainer->Add(fhConvAsymMCConversion) ;
506
507 fhConvPtMCConversion = new TH2F
508 ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
509 fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
510 fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
511 outputContainer->Add(fhConvPtMCConversion) ;
512
513 fhConvDispersionMCConversion = new TH2F
514 ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
515 fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
516 fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
517 outputContainer->Add(fhConvDispersionMCConversion) ;
518
519 fhConvM02MCConversion = new TH2F
520 ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
521 fhConvM02MCConversion->SetYTitle("M02 cluster 1");
522 fhConvM02MCConversion->SetXTitle("M02 cluster 2");
523 outputContainer->Add(fhConvM02MCConversion) ;
524
525 fhConvDeltaEtaMCAntiNeutron = new TH2F
526 ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
527 fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
528 fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
529 outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
530
531 fhConvDeltaPhiMCAntiNeutron = new TH2F
532 ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
533 fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
534 fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
535 outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
536
537 fhConvDeltaEtaPhiMCAntiNeutron = new TH2F
538 ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
539 fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
540 fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
541 outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
542
543 fhConvAsymMCAntiNeutron = new TH2F
544 ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
545 fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
546 fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
547 outputContainer->Add(fhConvAsymMCAntiNeutron) ;
548
549 fhConvPtMCAntiNeutron = new TH2F
550 ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
551 fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
552 fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
553 outputContainer->Add(fhConvPtMCAntiNeutron) ;
554
555 fhConvDispersionMCAntiNeutron = new TH2F
556 ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
557 fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
558 fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
559 outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
560
561 fhConvM02MCAntiNeutron = new TH2F
562 ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
563 fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
564 fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
565 outputContainer->Add(fhConvM02MCAntiNeutron) ;
566
567 fhConvDeltaEtaMCAntiProton = new TH2F
568 ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
569 fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
570 fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
571 outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
572
573 fhConvDeltaPhiMCAntiProton = new TH2F
574 ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
575 fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
576 fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
577 outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
578
579 fhConvDeltaEtaPhiMCAntiProton = new TH2F
580 ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
581 fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
582 fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
583 outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
584
585 fhConvAsymMCAntiProton = new TH2F
586 ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
587 fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
588 fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
589 outputContainer->Add(fhConvAsymMCAntiProton) ;
590
591 fhConvPtMCAntiProton = new TH2F
592 ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
593 fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
594 fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
595 outputContainer->Add(fhConvPtMCAntiProton) ;
596
597 fhConvDispersionMCAntiProton = new TH2F
598 ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
599 fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
600 fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
601 outputContainer->Add(fhConvDispersionMCAntiProton) ;
602
603 fhConvM02MCAntiProton = new TH2F
604 ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
605 fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
606 fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
607 outputContainer->Add(fhConvM02MCAntiProton) ;
608
609 fhConvDeltaEtaMCString = new TH2F
610 ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
611 fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
612 fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
613 outputContainer->Add(fhConvDeltaEtaMCString) ;
614
615 fhConvDeltaPhiMCString = new TH2F
616 ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
617 fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
618 fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
619 outputContainer->Add(fhConvDeltaPhiMCString) ;
620
621 fhConvDeltaEtaPhiMCString = new TH2F
622 ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
623 fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
624 fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
625 outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
626
627 fhConvAsymMCString = new TH2F
628 ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
629 fhConvAsymMCString->SetYTitle("Asymmetry");
630 fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
631 outputContainer->Add(fhConvAsymMCString) ;
632
633 fhConvPtMCString = new TH2F
634 ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
635 fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
636 fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
637 outputContainer->Add(fhConvPtMCString) ;
638
639 fhConvDispersionMCString = new TH2F
640 ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
641 fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
642 fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
643 outputContainer->Add(fhConvDispersionMCString) ;
644
645 fhConvM02MCString = new TH2F
646 ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
647 fhConvM02MCString->SetYTitle("M02 cluster 1");
648 fhConvM02MCString->SetXTitle("M02 cluster 2");
649 outputContainer->Add(fhConvM02MCString) ;
650 }
6175da48 651
477d6cee 652 }//Histos with MC
0c1383b5 653
477d6cee 654 return outputContainer ;
655
1c5acb87 656}
657
6639984f 658//____________________________________________________________________________
659void AliAnaPhoton::Init()
660{
661
662 //Init
663 //Do some checks
1e86c71e 664 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 665 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 666 abort();
667 }
1e86c71e 668 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 669 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 670 abort();
671 }
672
673}
674
675
1c5acb87 676//____________________________________________________________________________
677void AliAnaPhoton::InitParameters()
678{
679
680 //Initialize the parameters of the analysis.
a3aebfff 681 AddToHistogramsName("AnaPhoton_");
682
6175da48 683 fCalorimeter = "EMCAL" ;
684 fMinDist = 2.;
685 fMinDist2 = 4.;
686 fMinDist3 = 5.;
687 fMassCut = 0.03; //30 MeV
1e86c71e 688
4cf55759 689 fTimeCutMin = -1;
690 fTimeCutMax = 9999999;
6175da48 691 fNCellsCut = 0;
2ac125bf 692
1e86c71e 693 fRejectTrackMatch = kTRUE ;
694 fCheckConversion = kFALSE;
20218aea 695 fRemoveConvertedPair = kFALSE;
1e86c71e 696 fAddConvertedPairsToAOD = kFALSE;
697
1c5acb87 698}
699
700//__________________________________________________________________
701void AliAnaPhoton::MakeAnalysisFillAOD()
702{
f8006433 703 //Do photon analysis and fill aods
f37fa8d2 704
6175da48 705 //Get the vertex
5025c139 706 Double_t v[3] = {0,0,0}; //vertex ;
707 GetReader()->GetVertex(v);
f8006433 708
f37fa8d2 709 //Select the Calorimeter of the photon
c8fe2783 710 TObjArray * pl = 0x0;
477d6cee 711 if(fCalorimeter == "PHOS")
be518ab0 712 pl = GetPHOSClusters();
477d6cee 713 else if (fCalorimeter == "EMCAL")
be518ab0 714 pl = GetEMCALClusters();
5ae09196 715
716 if(!pl) {
717 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
718 return;
719 }
0ae57829 720
6175da48 721 //Init arrays, variables, get number of clusters
1e86c71e 722 TLorentzVector mom, mom2 ;
723 Int_t nCaloClusters = pl->GetEntriesFast();
6175da48 724 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
20218aea 725 Bool_t * indexConverted = 0x0;
726 if(fCheckConversion){
727 indexConverted = new Bool_t[nCaloClusters];
728 for (Int_t i = 0; i < nCaloClusters; i++)
729 indexConverted[i] = kFALSE;
730 }
731
6175da48 732 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
733
734 //----------------------------------------------------
735 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
736 //----------------------------------------------------
737 // Loop on clusters
1e86c71e 738 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
739
0ae57829 740 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
741 //printf("calo %d, %f\n",icalo,calo->E());
f8006433 742
743 //Get the index where the cluster comes, to retrieve the corresponding vertex
c8fe2783 744 Int_t evtIndex = 0 ;
745 if (GetMixedEvent()) {
746 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
5025c139 747 //Get the vertex and check it is not too large in z
96539743 748 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
c8fe2783 749 }
f8006433 750
f37fa8d2 751 //Cluster selection, not charged, with photon id and in fiducial cut
233e0df8 752
f37fa8d2 753 //Input from second AOD?
f8006433 754 //Int_t input = 0;
be518ab0 755 // if (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo)
f37fa8d2 756 // input = 1 ;
be518ab0 757 // else if(fCalorimeter == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= icalo)
f37fa8d2 758 // input = 1;
233e0df8 759
f37fa8d2 760 //Get Momentum vector,
f8006433 761 //if (input == 0)
762 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
763 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
764 else{
765 Double_t vertex[]={0,0,0};
766 calo->GetMomentum(mom,vertex) ;
767 }
f8006433 768
f37fa8d2 769 // else if(input == 1)
770 // calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
c8fe2783 771
6175da48 772 //--------------------------------------
773 // Cluster selection
774 //--------------------------------------
775 if(GetDebug() > 2)
1510eee3 776 printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
777 GetReader()->GetEventNumber(),
778 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
6175da48 779
780 //.......................................
f37fa8d2 781 //If too small or big pt, skip it
ba1eeb1f 782 if(mom.E() < GetMinPt() || mom.E() > GetMaxPt() ) continue ;
783 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",icalo);
c8fe2783 784
6175da48 785 //.......................................
786 // TOF cut, BE CAREFUL WITH THIS CUT
787 Double_t tof = calo->GetTOF()*1e9;
c8fe2783 788 if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
6175da48 789 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",icalo);
790
791 //.......................................
0ae57829 792 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
6175da48 793 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",icalo);
c8fe2783 794
6175da48 795 //.......................................
f37fa8d2 796 //Check acceptance selection
ff45398a 797 if(IsFiducialCutOn()){
798 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 799 if(! in ) continue ;
800 }
6175da48 801 if(GetDebug() > 2) printf("Fiducial cut passed \n");
c8fe2783 802
6175da48 803 //.......................................
804 //Skip matched clusters with tracks
805 if(fRejectTrackMatch){
806 if(IsTrackMatched(calo)) {
807 if(GetDebug() > 2) printf("\t Reject matched clusters\n");
808 continue ;
809 }
810 else
811 if(GetDebug() > 2) printf(" matching cut passed cut passed \n");
812 }// reject matched clusters
477d6cee 813
6175da48 814 //.......................................
f37fa8d2 815 //Check Distance to Bad channel, set bit.
c8fe2783 816 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 817 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
6175da48 818 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
477d6cee 819 continue ;
6175da48 820 }
821 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
477d6cee 822
6175da48 823 if(GetDebug() > 0)
1510eee3 824 printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
825 GetReader()->GetEventNumber(),
6175da48 826 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
477d6cee 827
6175da48 828
829 //----------------------------
830 //Create AOD for analysis
831 //----------------------------
832 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
833
834 //...............................................
835 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
836 Int_t label = calo->GetLabel();
837 aodph.SetLabel(label);
838 //aodph.SetInputFileIndex(input);
839 aodph.SetCaloLabel(calo->GetID(),-1);
840 aodph.SetDetector(fCalorimeter);
841 //printf("Index %d, Id %d\n",icalo, calo->GetID());
842
843 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
844
845 //...............................................
846 //Set bad channel distance bit
f37fa8d2 847 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
477d6cee 848 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
f37fa8d2 849 else aodph.SetDistToBad(0) ;
af7b3903 850 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
c8fe2783 851
6175da48 852 //...............................................
853 //Set number of cells in this cluster
854 //Temporary patch FIXME
855 aodph.SetBtag(calo->GetNCells());
856 // MEFIX
857
858 //-------------------------------------
f37fa8d2 859 //PID selection or bit setting
6175da48 860 //-------------------------------------
861 // MC
477d6cee 862 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
f37fa8d2 863 //Get most probable PID, check PID weights (in MC this option is mandatory)
c8fe2783 864 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
a3aebfff 865 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
f37fa8d2 866 //If primary is not photon, skip it.
477d6cee 867 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
6175da48 868 }
869 //...............................................
870 // Data, PID check on
477d6cee 871 else if(IsCaloPIDOn()){
f37fa8d2 872 //Get most probable PID, 2 options check PID weights
873 //or redo PID, recommended option for EMCal.
477d6cee 874 if(!IsCaloPIDRecalculationOn())
c8fe2783 875 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 876 else
c8fe2783 877 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
477d6cee 878
a3aebfff 879 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 880
f37fa8d2 881 //If cluster does not pass pid, not photon, skip it.
477d6cee 882 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
883
884 }
6175da48 885 //...............................................
886 // Data, PID check off
477d6cee 887 else{
f37fa8d2 888 //Set PID bits for later selection (AliAnaPi0 for example)
889 //GetPDG already called in SetPIDBits.
f2ccb5b8 890 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
a3aebfff 891 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 892 }
893
a3aebfff 894 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
477d6cee 895
6175da48 896 //--------------------------------------------------------------------------------------
f37fa8d2 897 //Play with the MC stack if available
6175da48 898 //--------------------------------------------------------------------------------------
899
f37fa8d2 900 //Check origin of the candidates
477d6cee 901 if(IsDataMC()){
c8fe2783 902 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
591cc579 903 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 904 }//Work with stack also
905
6175da48 906 //--------------------------------------------------------------------------------------
907 // Conversions pairs analysis
f37fa8d2 908 // Check if cluster comes from a conversion in the material in front of the calorimeter
909 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
6175da48 910 //--------------------------------------------------------------------------------------
911
912 // Do analysis only if there are more than one cluster
20218aea 913 if( nCaloClusters > 1 && fCheckConversion){
c8fe2783 914 Bool_t bConverted = kFALSE;
915 Int_t id2 = -1;
1e86c71e 916
f37fa8d2 917 //Check if set previously as converted couple, if so skip its use.
20218aea 918 if (indexConverted[icalo]) continue;
1e86c71e 919
6175da48 920 // Second cluster loop
c8fe2783 921 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
f37fa8d2 922 //Check if set previously as converted couple, if so skip its use.
c8fe2783 923 if (indexConverted[jcalo]) continue;
f37fa8d2 924 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
0ae57829 925 AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
6175da48 926
927 //Mixed event, get index of event
c8fe2783 928 Int_t evtIndex2 = 0 ;
929 if (GetMixedEvent()) {
930 evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
f8006433 931
6175da48 932 }
933
934 //Get kinematics of second cluster
f8006433 935 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
8cdc266d 936 calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
f8006433 937 else{
938 Double_t vertex[]={0,0,0};
8cdc266d 939 calo2->GetMomentum(mom2,vertex) ;
f8006433 940 }
941
f37fa8d2 942 //Check only certain regions
c8fe2783 943 Bool_t in2 = kTRUE;
944 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
945 if(!in2) continue;
946
6175da48 947 //................................................
f37fa8d2 948 //Get mass of pair, if small, take this pair.
41121cfe 949 Float_t pairM = (mom+mom2).M();
950 //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
951 if(pairM < fMassCut){
952 aodph.SetTagged(kFALSE);
c8fe2783 953 id2 = calo2->GetID();
6175da48 954 indexConverted[icalo]=kTRUE;
c8fe2783 955 indexConverted[jcalo]=kTRUE;
6175da48 956
41121cfe 957 Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
958 Float_t dPhi = mom.Phi()-mom2.Phi();
959 Float_t dEta = mom.Eta()-mom2.Eta();
960
8cdc266d 961 if(GetDebug() > 2)
41121cfe 962 printf("AliAnaPhoton::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",
963 pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
8cdc266d 964 calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
965 id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
6175da48 966
967 //...............................................
968 //Fill few histograms with kinematics of the pair
969 //FIXME, move all this to MakeAnalysisFillHistograms ...
41121cfe 970
971 fhConvDeltaEta ->Fill( pairM, dPhi );
972 fhConvDeltaPhi ->Fill( pairM, dEta );
973 fhConvAsym ->Fill( pairM, asymmetry );
974 fhConvDeltaEtaPhi->Fill( dEta , dPhi );
975 fhConvPt ->Fill( pairM, (mom+mom2).Pt());
6175da48 976
977 //...............................................
978 //Select pairs in a eta-phi window
41121cfe 979 if(TMath::Abs(dEta) < fConvDEtaCut &&
980 TMath::Abs(dPhi) < fConvDPhiMaxCut &&
981 TMath::Abs(dPhi) > fConvDPhiMinCut &&
982 asymmetry < fConvAsymCut ){
983 bConverted = kTRUE;
984 }
985 //printf("Accepted? %d\n",bConverted);
6175da48 986 //...........................................
987 //Fill more histograms, simulated data
988 //FIXME, move all this to MakeAnalysisFillHistograms ...
989 if(IsDataMC()){
990
991 //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
992 Int_t ancPDG = 0;
993 Int_t ancStatus = 0;
994 TLorentzVector momentum;
995 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
996 GetReader(), ancPDG, ancStatus, momentum);
997
998 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
999 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1000
1001 Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
1002 if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
1003 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
41121cfe 1004 fhConvDeltaEtaMCConversion ->Fill( pairM, dEta );
1005 fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi );
1006 fhConvAsymMCConversion ->Fill( pairM, asymmetry );
1007 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi );
1008 fhConvPtMCConversion ->Fill( pairM, (mom+mom2).Pt());
6175da48 1009 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1010 fhConvM02MCConversion ->Fill( calo->GetM02(), calo2->GetM02());
1011
1012 }
1013 }
1014 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
1015 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
41121cfe 1016 fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta );
1017 fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi );
1018 fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry );
1019 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi );
1020 fhConvPtMCAntiNeutron ->Fill( pairM, (mom+mom2).Pt());
6175da48 1021 fhConvDispersionMCAntiNeutron ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1022 fhConvM02MCAntiNeutron ->Fill( calo->GetM02(), calo2->GetM02());
1023 }
1024 }
1025 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
1026 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
41121cfe 1027 fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta );
1028 fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi );
1029 fhConvAsymMCAntiProton ->Fill( pairM, asymmetry );
1030 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi );
1031 fhConvPtMCAntiProton ->Fill( pairM, (mom+mom2).Pt());
6175da48 1032 fhConvDispersionMCAntiProton ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1033 fhConvM02MCAntiProton ->Fill( calo->GetM02(), calo2->GetM02());
1034 }
1035 }
1036
1037 //Pairs coming from fragmenting pairs.
1038 if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
41121cfe 1039 fhConvDeltaEtaMCString ->Fill( pairM, dPhi);
1040 fhConvDeltaPhiMCString ->Fill( pairM, dPhi);
1041 fhConvAsymMCString ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
1042 fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
1043 fhConvPtMCString ->Fill( pairM, (mom+mom2).Pt());
6175da48 1044 fhConvDispersionMCString ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1045 fhConvM02MCString ->Fill( calo->GetM02(), calo2->GetM02());
1046 }
1047
1048 }// Data MC
1049
c8fe2783 1050 break;
1051 }
1e86c71e 1052
c8fe2783 1053 }//Mass loop
1e86c71e 1054
6175da48 1055 //..........................................................................................................
1056 //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
c8fe2783 1057 if(bConverted){
6175da48 1058 //Add to AOD
c8fe2783 1059 if(fAddConvertedPairsToAOD){
f37fa8d2 1060 //Create AOD of pair analysis
c8fe2783 1061 TLorentzVector mpair = mom+mom2;
1062 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
1063 aodpair.SetLabel(aodph.GetLabel());
f8006433 1064 //aodpair.SetInputFileIndex(input);
c8fe2783 1065
f37fa8d2 1066 //printf("Index %d, Id %d\n",icalo, calo->GetID());
1067 //Set the indeces of the original caloclusters
c8fe2783 1068 aodpair.SetCaloLabel(calo->GetID(),id2);
1069 aodpair.SetDetector(fCalorimeter);
1070 aodpair.SetPdg(aodph.GetPdg());
1071 aodpair.SetTag(aodph.GetTag());
6175da48 1072 aodpair.SetTagged(kTRUE);
f37fa8d2 1073 //Add AOD with pair object to aod branch
c8fe2783 1074 AddAODParticle(aodpair);
f37fa8d2 1075 //printf("\t \t both added pair\n");
c8fe2783 1076 }
1077
f37fa8d2 1078 //Do not add the current calocluster
20218aea 1079 if(fRemoveConvertedPair) continue;
6175da48 1080 else {
1081 //printf("TAGGED\n");
1082 //Tag this cluster as likely conversion
1083 aodph.SetTagged(kTRUE);
1084 }
c8fe2783 1085 }//converted pair
1086 }//check conversion
f37fa8d2 1087 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 1088
6175da48 1089 //FIXME, this to MakeAnalysisFillHistograms ...
1090 fhNCellsPt->Fill(aodph.Pt(),calo->GetNCells());
1091
f37fa8d2 1092 //Add AOD with photon object to aod branch
477d6cee 1093 AddAODParticle(aodph);
1094
1095 }//loop
1096
4a745797 1097 delete [] indexConverted;
1098
f37fa8d2 1099 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 1100
1c5acb87 1101}
1102
1103//__________________________________________________________________
1104void AliAnaPhoton::MakeAnalysisFillHistograms()
1105{
6175da48 1106 //Fill histograms
f8006433 1107
6175da48 1108 //-------------------------------------------------------------------
591cc579 1109 // Access MC information in stack if requested, check that it exists.
1110 AliStack * stack = 0x0;
1111 TParticle * primary = 0x0;
1112 TClonesArray * mcparticles0 = 0x0;
f8006433 1113 //TClonesArray * mcparticles1 = 0x0;
591cc579 1114 AliAODMCParticle * aodprimary = 0x0;
1115 if(IsDataMC()){
1116
1117 if(GetReader()->ReadStack()){
1118 stack = GetMCStack() ;
1119 if(!stack) {
1120 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1121 abort();
1122 }
f8006433 1123
591cc579 1124 }
1125 else if(GetReader()->ReadAODMCParticles()){
f8006433 1126
591cc579 1127 //Get the list of MC particles
1128 mcparticles0 = GetReader()->GetAODMCParticles(0);
1129 if(!mcparticles0 && GetDebug() > 0) {
1130 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1131 }
f8006433 1132 // if(GetReader()->GetSecondInputAODTree()){
1133 // mcparticles1 = GetReader()->GetAODMCParticles(1);
1134 // if(!mcparticles1 && GetDebug() > 0) {
1135 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
1136 // }
1137 // }
591cc579 1138
1139 }
1140 }// is data and MC
1141
6175da48 1142
1143 // Get vertex
2244659d 1144 Double_t v[3] = {0,0,0}; //vertex ;
1145 GetReader()->GetVertex(v);
6175da48 1146 //fhVertex->Fill(v[0],v[1],v[2]);
1147 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1148
1149 //----------------------------------
1150 //Loop on stored AOD photons
123fc3bd 1151 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2244659d 1152 fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
123fc3bd 1153 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1154
1155 for(Int_t iaod = 0; iaod < naod ; iaod++){
c8fe2783 1156 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1157 Int_t pdg = ph->GetPdg();
1158
1159 if(GetDebug() > 3)
1160 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1161
1162 //If PID used, fill histos with photons in Calorimeter fCalorimeter
1163 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
1164 if(ph->GetDetector() != fCalorimeter) continue;
1165
1166 if(GetDebug() > 2)
1167 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1168
6175da48 1169 //................................
c8fe2783 1170 //Fill photon histograms
1171 Float_t ptcluster = ph->Pt();
1172 Float_t phicluster = ph->Phi();
1173 Float_t etacluster = ph->Eta();
1174 Float_t ecluster = ph->E();
1175
20218aea 1176 fhEPhoton ->Fill(ecluster);
c8fe2783 1177 fhPtPhoton ->Fill(ptcluster);
1178 fhPhiPhoton ->Fill(ptcluster,phicluster);
20218aea 1179 fhEtaPhoton ->Fill(ptcluster,etacluster);
1180 if(ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
1181 else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
1182
1183 if(fCheckConversion &&ph->IsTagged()){
1184 fhPtPhotonConv->Fill(ptcluster);
1185 if(ecluster > 0.5) fhEtaPhiPhotonConv ->Fill(etacluster, phicluster);
1186 else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
6175da48 1187 }
20218aea 1188
6175da48 1189 //.......................................
c8fe2783 1190 //Play with the MC data if available
1191 if(IsDataMC()){
1192
1193 Int_t tag =ph->GetTag();
1194
1195 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
f8006433 1196 {
1197 fhPtMCPhoton ->Fill(ptcluster);
1198 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
1199 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
1200
1201 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
1202 {
1203 fhPtConversion ->Fill(ptcluster);
1204 fhPhiConversion ->Fill(ptcluster,phicluster);
1205 fhEtaConversion ->Fill(ptcluster,etacluster);
6175da48 1206 if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
1207 if(ptcluster > 0.5)fhEtaPhiConversion ->Fill(etacluster,phicluster);
1208 else fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
f8006433 1209 }
1210
1211 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
1212 fhPtPrompt ->Fill(ptcluster);
1213 fhPhiPrompt ->Fill(ptcluster,phicluster);
1214 fhEtaPrompt ->Fill(ptcluster,etacluster);
1215 }
1216 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
1217 {
1218 fhPtFragmentation ->Fill(ptcluster);
1219 fhPhiFragmentation ->Fill(ptcluster,phicluster);
1220 fhEtaFragmentation ->Fill(ptcluster,etacluster);
1221 }
1222 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
1223 {
1224 fhPtISR ->Fill(ptcluster);
1225 fhPhiISR ->Fill(ptcluster,phicluster);
1226 fhEtaISR ->Fill(ptcluster,etacluster);
1227 }
1228 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
1229 {
1230 fhPtPi0Decay ->Fill(ptcluster);
1231 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
1232 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
1233 }
1234 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
1235 {
1236 fhPtOtherDecay ->Fill(ptcluster);
1237 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
1238 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
1239 }
1240 }
6175da48 1241 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
1242 {
1243 fhPtAntiNeutron ->Fill(ptcluster);
1244 fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
1245 fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
20218aea 1246 if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
41121cfe 1247
6175da48 1248 }
1249 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
1250 {
1251 fhPtAntiProton ->Fill(ptcluster);
1252 fhPhiAntiProton ->Fill(ptcluster,phicluster);
1253 fhEtaAntiProton ->Fill(ptcluster,etacluster);
20218aea 1254 if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
41121cfe 1255
6175da48 1256 }
c8fe2783 1257 else{
1258 fhPtUnknown ->Fill(ptcluster);
1259 fhPhiUnknown ->Fill(ptcluster,phicluster);
1260 fhEtaUnknown ->Fill(ptcluster,etacluster);
20218aea 1261 if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
41121cfe 1262
c8fe2783 1263
f8006433 1264 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1265 // ph->GetLabel(),ph->Pt());
1266 // for(Int_t i = 0; i < 20; i++) {
1267 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1268 // }
1269 // printf("\n");
1270
c8fe2783 1271 }
1272
6175da48 1273 //....................................................................
c8fe2783 1274 // Access MC information in stack if requested, check that it exists.
1275 Int_t label =ph->GetLabel();
1276 if(label < 0) {
1277 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1278 continue;
1279 }
1280
1281 Float_t eprim = 0;
1282 Float_t ptprim = 0;
1283 if(GetReader()->ReadStack()){
1284
1285 if(label >= stack->GetNtrack()) {
f8006433 1286 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1287 continue ;
c8fe2783 1288 }
1289
1290 primary = stack->Particle(label);
1291 if(!primary){
f8006433 1292 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1293 continue;
c8fe2783 1294 }
1295 eprim = primary->Energy();
1296 ptprim = primary->Pt();
1297
1298 }
1299 else if(GetReader()->ReadAODMCParticles()){
1300 //Check which is the input
1301 if(ph->GetInputFileIndex() == 0){
f8006433 1302 if(!mcparticles0) continue;
1303 if(label >= mcparticles0->GetEntriesFast()) {
1304 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1305 label, mcparticles0->GetEntriesFast());
1306 continue ;
1307 }
1308 //Get the particle
1309 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1310
c8fe2783 1311 }
f8006433 1312// else {//Second input
1313// if(!mcparticles1) continue;
1314// if(label >= mcparticles1->GetEntriesFast()) {
1315// if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1316// label, mcparticles1->GetEntriesFast());
1317// continue ;
1318// }
1319// //Get the particle
1320// aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
1321//
1322// }//second input
c8fe2783 1323
1324 if(!aodprimary){
f8006433 1325 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1326 continue;
c8fe2783 1327 }
1328
1329 eprim = aodprimary->E();
1330 ptprim = aodprimary->Pt();
1331
1332 }
1333
1334 fh2E ->Fill(ecluster, eprim);
1335 fh2Pt ->Fill(ptcluster, ptprim);
1336 fhDeltaE ->Fill(eprim-ecluster);
1337 fhDeltaPt->Fill(ptprim-ptcluster);
1338 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
1339 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
1340
1341 }//Histograms with MC
1342
1343 }// aod loop
591cc579 1344
1c5acb87 1345}
1346
1347
1348//__________________________________________________________________
1349void AliAnaPhoton::Print(const Option_t * opt) const
1350{
477d6cee 1351 //Print some relevant parameters set for the analysis
1352
1353 if(! opt)
1354 return;
1355
1356 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1357 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 1358
477d6cee 1359 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1360 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1361 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1362 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 1363 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 1364 printf("Check Pair Conversion = %d\n",fCheckConversion);
1365 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
1366 printf("Conversion pair mass cut = %f\n",fMassCut);
41121cfe 1367 printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
1368 fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
4cf55759 1369 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 1370 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 1371 printf(" \n") ;
1c5acb87 1372
1373}