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