Fixed a bug which caused a wrong filling of histograms
[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
fb516de2 161 fhNCellsPt = new TH2F ("hNCellsPt","# of cells in cluster vs E of clusters", nptbins,ptmin, ptmax, 100,0,100);
ba1eeb1f 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
658//____________________________________________________________________________
6639984f 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
676//____________________________________________________________________________
1c5acb87 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()));
fb516de2 903 if(GetDebug() > 0)
904 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 905 }//Work with stack also
906
6175da48 907 //--------------------------------------------------------------------------------------
908 // Conversions pairs analysis
f37fa8d2 909 // Check if cluster comes from a conversion in the material in front of the calorimeter
910 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
6175da48 911 //--------------------------------------------------------------------------------------
912
913 // Do analysis only if there are more than one cluster
20218aea 914 if( nCaloClusters > 1 && fCheckConversion){
c8fe2783 915 Bool_t bConverted = kFALSE;
916 Int_t id2 = -1;
1e86c71e 917
f37fa8d2 918 //Check if set previously as converted couple, if so skip its use.
20218aea 919 if (indexConverted[icalo]) continue;
1e86c71e 920
6175da48 921 // Second cluster loop
c8fe2783 922 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
f37fa8d2 923 //Check if set previously as converted couple, if so skip its use.
c8fe2783 924 if (indexConverted[jcalo]) continue;
f37fa8d2 925 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
0ae57829 926 AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
6175da48 927
928 //Mixed event, get index of event
c8fe2783 929 Int_t evtIndex2 = 0 ;
930 if (GetMixedEvent()) {
931 evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
f8006433 932
6175da48 933 }
934
935 //Get kinematics of second cluster
f8006433 936 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
8cdc266d 937 calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
f8006433 938 else{
939 Double_t vertex[]={0,0,0};
8cdc266d 940 calo2->GetMomentum(mom2,vertex) ;
f8006433 941 }
942
f37fa8d2 943 //Check only certain regions
c8fe2783 944 Bool_t in2 = kTRUE;
945 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
946 if(!in2) continue;
947
6175da48 948 //................................................
f37fa8d2 949 //Get mass of pair, if small, take this pair.
41121cfe 950 Float_t pairM = (mom+mom2).M();
951 //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
952 if(pairM < fMassCut){
953 aodph.SetTagged(kFALSE);
c8fe2783 954 id2 = calo2->GetID();
6175da48 955 indexConverted[icalo]=kTRUE;
c8fe2783 956 indexConverted[jcalo]=kTRUE;
6175da48 957
41121cfe 958 Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
959 Float_t dPhi = mom.Phi()-mom2.Phi();
960 Float_t dEta = mom.Eta()-mom2.Eta();
961
8cdc266d 962 if(GetDebug() > 2)
41121cfe 963 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",
964 pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
8cdc266d 965 calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
966 id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
6175da48 967
968 //...............................................
969 //Fill few histograms with kinematics of the pair
970 //FIXME, move all this to MakeAnalysisFillHistograms ...
41121cfe 971
972 fhConvDeltaEta ->Fill( pairM, dPhi );
973 fhConvDeltaPhi ->Fill( pairM, dEta );
974 fhConvAsym ->Fill( pairM, asymmetry );
975 fhConvDeltaEtaPhi->Fill( dEta , dPhi );
976 fhConvPt ->Fill( pairM, (mom+mom2).Pt());
6175da48 977
978 //...............................................
979 //Select pairs in a eta-phi window
41121cfe 980 if(TMath::Abs(dEta) < fConvDEtaCut &&
981 TMath::Abs(dPhi) < fConvDPhiMaxCut &&
982 TMath::Abs(dPhi) > fConvDPhiMinCut &&
983 asymmetry < fConvAsymCut ){
984 bConverted = kTRUE;
985 }
986 //printf("Accepted? %d\n",bConverted);
6175da48 987 //...........................................
988 //Fill more histograms, simulated data
989 //FIXME, move all this to MakeAnalysisFillHistograms ...
990 if(IsDataMC()){
991
992 //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
993 Int_t ancPDG = 0;
994 Int_t ancStatus = 0;
995 TLorentzVector momentum;
996 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
997 GetReader(), ancPDG, ancStatus, momentum);
998
999 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1000 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1001
1002 Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
1003 if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
1004 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
41121cfe 1005 fhConvDeltaEtaMCConversion ->Fill( pairM, dEta );
1006 fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi );
1007 fhConvAsymMCConversion ->Fill( pairM, asymmetry );
1008 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi );
1009 fhConvPtMCConversion ->Fill( pairM, (mom+mom2).Pt());
6175da48 1010 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1011 fhConvM02MCConversion ->Fill( calo->GetM02(), calo2->GetM02());
1012
1013 }
1014 }
1015 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
1016 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
41121cfe 1017 fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta );
1018 fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi );
1019 fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry );
1020 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi );
1021 fhConvPtMCAntiNeutron ->Fill( pairM, (mom+mom2).Pt());
6175da48 1022 fhConvDispersionMCAntiNeutron ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1023 fhConvM02MCAntiNeutron ->Fill( calo->GetM02(), calo2->GetM02());
1024 }
1025 }
1026 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
1027 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
41121cfe 1028 fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta );
1029 fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi );
1030 fhConvAsymMCAntiProton ->Fill( pairM, asymmetry );
1031 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi );
1032 fhConvPtMCAntiProton ->Fill( pairM, (mom+mom2).Pt());
6175da48 1033 fhConvDispersionMCAntiProton ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1034 fhConvM02MCAntiProton ->Fill( calo->GetM02(), calo2->GetM02());
1035 }
1036 }
1037
1038 //Pairs coming from fragmenting pairs.
1039 if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
41121cfe 1040 fhConvDeltaEtaMCString ->Fill( pairM, dPhi);
1041 fhConvDeltaPhiMCString ->Fill( pairM, dPhi);
1042 fhConvAsymMCString ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
1043 fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
1044 fhConvPtMCString ->Fill( pairM, (mom+mom2).Pt());
6175da48 1045 fhConvDispersionMCString ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1046 fhConvM02MCString ->Fill( calo->GetM02(), calo2->GetM02());
1047 }
1048
1049 }// Data MC
1050
c8fe2783 1051 break;
1052 }
1e86c71e 1053
c8fe2783 1054 }//Mass loop
1e86c71e 1055
6175da48 1056 //..........................................................................................................
1057 //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
c8fe2783 1058 if(bConverted){
6175da48 1059 //Add to AOD
c8fe2783 1060 if(fAddConvertedPairsToAOD){
f37fa8d2 1061 //Create AOD of pair analysis
c8fe2783 1062 TLorentzVector mpair = mom+mom2;
1063 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
1064 aodpair.SetLabel(aodph.GetLabel());
f8006433 1065 //aodpair.SetInputFileIndex(input);
c8fe2783 1066
f37fa8d2 1067 //printf("Index %d, Id %d\n",icalo, calo->GetID());
1068 //Set the indeces of the original caloclusters
c8fe2783 1069 aodpair.SetCaloLabel(calo->GetID(),id2);
1070 aodpair.SetDetector(fCalorimeter);
1071 aodpair.SetPdg(aodph.GetPdg());
1072 aodpair.SetTag(aodph.GetTag());
6175da48 1073 aodpair.SetTagged(kTRUE);
f37fa8d2 1074 //Add AOD with pair object to aod branch
c8fe2783 1075 AddAODParticle(aodpair);
f37fa8d2 1076 //printf("\t \t both added pair\n");
c8fe2783 1077 }
1078
f37fa8d2 1079 //Do not add the current calocluster
20218aea 1080 if(fRemoveConvertedPair) continue;
6175da48 1081 else {
1082 //printf("TAGGED\n");
1083 //Tag this cluster as likely conversion
1084 aodph.SetTagged(kTRUE);
1085 }
c8fe2783 1086 }//converted pair
1087 }//check conversion
f37fa8d2 1088 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 1089
6175da48 1090 //FIXME, this to MakeAnalysisFillHistograms ...
1091 fhNCellsPt->Fill(aodph.Pt(),calo->GetNCells());
1092
f37fa8d2 1093 //Add AOD with photon object to aod branch
477d6cee 1094 AddAODParticle(aodph);
1095
1096 }//loop
1097
4a745797 1098 delete [] indexConverted;
1099
f37fa8d2 1100 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 1101
1c5acb87 1102}
1103
1104//__________________________________________________________________
1105void AliAnaPhoton::MakeAnalysisFillHistograms()
1106{
6175da48 1107 //Fill histograms
f8006433 1108
6175da48 1109 //-------------------------------------------------------------------
591cc579 1110 // Access MC information in stack if requested, check that it exists.
1111 AliStack * stack = 0x0;
1112 TParticle * primary = 0x0;
1113 TClonesArray * mcparticles0 = 0x0;
f8006433 1114 //TClonesArray * mcparticles1 = 0x0;
591cc579 1115 AliAODMCParticle * aodprimary = 0x0;
1116 if(IsDataMC()){
1117
1118 if(GetReader()->ReadStack()){
1119 stack = GetMCStack() ;
1120 if(!stack) {
1121 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1122 abort();
1123 }
f8006433 1124
591cc579 1125 }
1126 else if(GetReader()->ReadAODMCParticles()){
f8006433 1127
591cc579 1128 //Get the list of MC particles
1129 mcparticles0 = GetReader()->GetAODMCParticles(0);
1130 if(!mcparticles0 && GetDebug() > 0) {
1131 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1132 }
f8006433 1133 // if(GetReader()->GetSecondInputAODTree()){
1134 // mcparticles1 = GetReader()->GetAODMCParticles(1);
1135 // if(!mcparticles1 && GetDebug() > 0) {
1136 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
1137 // }
1138 // }
591cc579 1139
1140 }
1141 }// is data and MC
1142
6175da48 1143
1144 // Get vertex
2244659d 1145 Double_t v[3] = {0,0,0}; //vertex ;
1146 GetReader()->GetVertex(v);
6175da48 1147 //fhVertex->Fill(v[0],v[1],v[2]);
1148 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1149
1150 //----------------------------------
1151 //Loop on stored AOD photons
123fc3bd 1152 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2244659d 1153 fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
123fc3bd 1154 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1155
1156 for(Int_t iaod = 0; iaod < naod ; iaod++){
c8fe2783 1157 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1158 Int_t pdg = ph->GetPdg();
1159
1160 if(GetDebug() > 3)
1161 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1162
1163 //If PID used, fill histos with photons in Calorimeter fCalorimeter
1164 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
1165 if(ph->GetDetector() != fCalorimeter) continue;
1166
1167 if(GetDebug() > 2)
1168 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1169
6175da48 1170 //................................
c8fe2783 1171 //Fill photon histograms
1172 Float_t ptcluster = ph->Pt();
1173 Float_t phicluster = ph->Phi();
1174 Float_t etacluster = ph->Eta();
1175 Float_t ecluster = ph->E();
1176
20218aea 1177 fhEPhoton ->Fill(ecluster);
c8fe2783 1178 fhPtPhoton ->Fill(ptcluster);
1179 fhPhiPhoton ->Fill(ptcluster,phicluster);
20218aea 1180 fhEtaPhoton ->Fill(ptcluster,etacluster);
1181 if(ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
1182 else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
1183
1184 if(fCheckConversion &&ph->IsTagged()){
1185 fhPtPhotonConv->Fill(ptcluster);
1186 if(ecluster > 0.5) fhEtaPhiPhotonConv ->Fill(etacluster, phicluster);
1187 else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
6175da48 1188 }
20218aea 1189
6175da48 1190 //.......................................
c8fe2783 1191 //Play with the MC data if available
1192 if(IsDataMC()){
1193
1194 Int_t tag =ph->GetTag();
fb516de2 1195
c8fe2783 1196 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
f8006433 1197 {
1198 fhPtMCPhoton ->Fill(ptcluster);
1199 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
1200 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
1201
1202 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
1203 {
1204 fhPtConversion ->Fill(ptcluster);
1205 fhPhiConversion ->Fill(ptcluster,phicluster);
1206 fhEtaConversion ->Fill(ptcluster,etacluster);
6175da48 1207 if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
1208 if(ptcluster > 0.5)fhEtaPhiConversion ->Fill(etacluster,phicluster);
1209 else fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
f8006433 1210 }
1211
1212 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
1213 fhPtPrompt ->Fill(ptcluster);
1214 fhPhiPrompt ->Fill(ptcluster,phicluster);
1215 fhEtaPrompt ->Fill(ptcluster,etacluster);
1216 }
1217 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
1218 {
1219 fhPtFragmentation ->Fill(ptcluster);
1220 fhPhiFragmentation ->Fill(ptcluster,phicluster);
1221 fhEtaFragmentation ->Fill(ptcluster,etacluster);
1222 }
1223 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
1224 {
1225 fhPtISR ->Fill(ptcluster);
1226 fhPhiISR ->Fill(ptcluster,phicluster);
1227 fhEtaISR ->Fill(ptcluster,etacluster);
1228 }
1229 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
1230 {
1231 fhPtPi0Decay ->Fill(ptcluster);
1232 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
1233 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
1234 }
1235 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
1236 {
1237 fhPtOtherDecay ->Fill(ptcluster);
1238 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
1239 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
1240 }
1241 }
6175da48 1242 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
1243 {
1244 fhPtAntiNeutron ->Fill(ptcluster);
1245 fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
1246 fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
20218aea 1247 if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
41121cfe 1248
6175da48 1249 }
1250 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
1251 {
1252 fhPtAntiProton ->Fill(ptcluster);
1253 fhPhiAntiProton ->Fill(ptcluster,phicluster);
1254 fhEtaAntiProton ->Fill(ptcluster,etacluster);
20218aea 1255 if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
41121cfe 1256
6175da48 1257 }
c8fe2783 1258 else{
1259 fhPtUnknown ->Fill(ptcluster);
1260 fhPhiUnknown ->Fill(ptcluster,phicluster);
1261 fhEtaUnknown ->Fill(ptcluster,etacluster);
20218aea 1262 if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
41121cfe 1263
c8fe2783 1264
f8006433 1265 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1266 // ph->GetLabel(),ph->Pt());
1267 // for(Int_t i = 0; i < 20; i++) {
1268 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1269 // }
1270 // printf("\n");
1271
c8fe2783 1272 }
1273
6175da48 1274 //....................................................................
c8fe2783 1275 // Access MC information in stack if requested, check that it exists.
1276 Int_t label =ph->GetLabel();
1277 if(label < 0) {
1278 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1279 continue;
1280 }
1281
1282 Float_t eprim = 0;
1283 Float_t ptprim = 0;
1284 if(GetReader()->ReadStack()){
1285
1286 if(label >= stack->GetNtrack()) {
f8006433 1287 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1288 continue ;
c8fe2783 1289 }
1290
1291 primary = stack->Particle(label);
1292 if(!primary){
f8006433 1293 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1294 continue;
c8fe2783 1295 }
1296 eprim = primary->Energy();
1297 ptprim = primary->Pt();
1298
1299 }
1300 else if(GetReader()->ReadAODMCParticles()){
1301 //Check which is the input
1302 if(ph->GetInputFileIndex() == 0){
f8006433 1303 if(!mcparticles0) continue;
1304 if(label >= mcparticles0->GetEntriesFast()) {
1305 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1306 label, mcparticles0->GetEntriesFast());
1307 continue ;
1308 }
1309 //Get the particle
1310 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1311
c8fe2783 1312 }
f8006433 1313// else {//Second input
1314// if(!mcparticles1) continue;
1315// if(label >= mcparticles1->GetEntriesFast()) {
1316// if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1317// label, mcparticles1->GetEntriesFast());
1318// continue ;
1319// }
1320// //Get the particle
1321// aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
1322//
1323// }//second input
c8fe2783 1324
1325 if(!aodprimary){
f8006433 1326 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1327 continue;
c8fe2783 1328 }
1329
1330 eprim = aodprimary->E();
1331 ptprim = aodprimary->Pt();
1332
1333 }
1334
1335 fh2E ->Fill(ecluster, eprim);
1336 fh2Pt ->Fill(ptcluster, ptprim);
1337 fhDeltaE ->Fill(eprim-ecluster);
1338 fhDeltaPt->Fill(ptprim-ptcluster);
1339 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
1340 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
1341
1342 }//Histograms with MC
1343
1344 }// aod loop
591cc579 1345
1c5acb87 1346}
1347
1348
1349//__________________________________________________________________
1350void AliAnaPhoton::Print(const Option_t * opt) const
1351{
477d6cee 1352 //Print some relevant parameters set for the analysis
1353
1354 if(! opt)
1355 return;
1356
1357 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1358 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 1359
477d6cee 1360 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1361 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1362 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1363 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 1364 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 1365 printf("Check Pair Conversion = %d\n",fCheckConversion);
1366 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
1367 printf("Conversion pair mass cut = %f\n",fMassCut);
41121cfe 1368 printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
1369 fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
4cf55759 1370 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 1371 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 1372 printf(" \n") ;
1c5acb87 1373
1374}