]> 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
adding implementations to extract event type and time stamp from RawReader
[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 *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $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.
22//
23// -- Author: Gustavo Conesa (LNF-INFN)
24//////////////////////////////////////////////////////////////////////////////
25
26
27// --- ROOT system ---
28#include <TH2F.h>
477d6cee 29#include <TClonesArray.h>
0c1383b5 30#include <TObjString.h>
477d6cee 31//#include <Riostream.h>
123fc3bd 32#include "TParticle.h"
1c5acb87 33
34// --- Analysis system ---
35#include "AliAnaPhoton.h"
36#include "AliCaloTrackReader.h"
123fc3bd 37#include "AliStack.h"
1c5acb87 38#include "AliCaloPID.h"
6639984f 39#include "AliMCAnalysisUtils.h"
ff45398a 40#include "AliFiducialCut.h"
477d6cee 41#include "AliAODCaloCluster.h"
591cc579 42#include "AliAODMCParticle.h"
1c5acb87 43
44ClassImp(AliAnaPhoton)
45
46//____________________________________________________________________________
47 AliAnaPhoton::AliAnaPhoton() :
48 AliAnaPartCorrBaseClass(), fCalorimeter(""),
a3aebfff 49 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
1e86c71e 50 fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
2ac125bf 51 fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
a3aebfff 52 fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
7175a03a 53 //MC
123fc3bd 54 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
591cc579 55 fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
1c5acb87 56 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
57 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
6639984f 58 fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
1c5acb87 59 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
60 fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
61 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
62 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
63{
64 //default ctor
65
66 //Initialize parameters
67 InitParameters();
68
69}
78219bac 70/*
1c5acb87 71//____________________________________________________________________________
72AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
73 AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
1e86c71e 74 fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
75 fRejectTrackMatch(g.fRejectTrackMatch),
76 fCheckConversion(g.fCheckConversion),fAddConvertedPairsToAOD(g.fAddConvertedPairsToAOD),
2ac125bf 77 fMassCut(g.fMassCut), fTimeCutMin(g.fTimeCutMin), fTimeCutMax(g.fTimeCutMax), fNCellsCut(g.fNCellsCut),
1e86c71e 78 fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
1c5acb87 79 //MC
123fc3bd 80 fhDeltaE(g.fhDeltaE), fhDeltaPt(g.fhDeltaPt),
81 fhRatioE(g.fhRatioE), fhRatioPt(g.fhRatioPt),
82 fh2E(g.fh2E), fh2Pt(g.fh2Pt),
591cc579 83 fhPtMCPhoton(g.fhPtMCPhoton),fhPhiMCPhoton(g.fhPhiMCPhoton),fhEtaMCPhoton(g.fhEtaMCPhoton),
1c5acb87 84 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
85 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
6639984f 86 fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
1c5acb87 87 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
88 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
89 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
6639984f 90 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
91
1c5acb87 92{
93 // cpy ctor
94
95}
96
97//_________________________________________________________________________
98AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
99{
100 // assignment operator
101
102 if(&g == this) return *this;
103
104 fCalorimeter = g.fCalorimeter ;
105 fMinDist = g.fMinDist;
106 fMinDist2 = g.fMinDist2;
107 fMinDist3 = g.fMinDist3;
4cf55759 108 fMassCut = g.fMassCut;
1e86c71e 109
4cf55759 110 fTimeCutMin = g.fTimeCutMin;
111 fTimeCutMax = g.fTimeCutMax;
2ac125bf 112 fNCellsCut = g.fNCellsCut;
113
1e86c71e 114 fRejectTrackMatch = g.fRejectTrackMatch;
115 fCheckConversion = g.fCheckConversion;
116 fAddConvertedPairsToAOD = g.fAddConvertedPairsToAOD;
117
6639984f 118 fhPtPhoton = g.fhPtPhoton ;
1c5acb87 119 fhPhiPhoton = g.fhPhiPhoton ;
120 fhEtaPhoton = g.fhEtaPhoton ;
123fc3bd 121
122 fhDeltaE = g.fhDeltaE;
123 fhDeltaPt = g.fhDeltaPt;
124 fhRatioE = g.fhRatioE;
125 fhRatioPt = g.fhRatioPt;
126 fh2E = g.fh2E;
127 fh2Pt = g.fh2Pt;
1c5acb87 128
591cc579 129 fhPtMCPhoton = g.fhPtMCPhoton ;
130 fhPhiMCPhoton = g.fhPhiMCPhoton ;
131 fhEtaMCPhoton = g.fhEtaMCPhoton ;
6639984f 132 fhPtPrompt = g.fhPtPrompt;
1c5acb87 133 fhPhiPrompt = g.fhPhiPrompt;
134 fhEtaPrompt = g.fhEtaPrompt;
6639984f 135 fhPtFragmentation = g.fhPtFragmentation;
1c5acb87 136 fhPhiFragmentation = g.fhPhiFragmentation;
137 fhEtaFragmentation = g.fhEtaFragmentation;
6639984f 138 fhPtISR = g.fhPtISR;
139 fhPhiISR = g.fhPhiISR;
140 fhEtaISR = g.fhEtaISR;
141 fhPtPi0Decay = g.fhPtPi0Decay;
1c5acb87 142 fhPhiPi0Decay = g.fhPhiPi0Decay;
143 fhEtaPi0Decay = g.fhEtaPi0Decay;
6639984f 144 fhPtOtherDecay = g.fhPtOtherDecay;
1c5acb87 145 fhPhiOtherDecay = g.fhPhiOtherDecay;
146 fhEtaOtherDecay = g.fhEtaOtherDecay;
6639984f 147 fhPtConversion = g. fhPtConversion;
1c5acb87 148 fhPhiConversion = g.fhPhiConversion;
149 fhEtaConversion = g.fhEtaConversion;
6639984f 150 fhPtUnknown = g.fhPtUnknown;
1c5acb87 151 fhPhiUnknown = g.fhPhiUnknown;
152 fhEtaUnknown = g.fhEtaUnknown;
153
154 return *this;
155
156}
78219bac 157*/
1c5acb87 158//____________________________________________________________________________
159AliAnaPhoton::~AliAnaPhoton()
160{
161 //dtor
162
163}
164
0c1383b5 165//________________________________________________________________________
166TObjString * AliAnaPhoton::GetAnalysisCuts()
167{
168 //Save parameters used for analysis
169 TString parList ; //this will be list of parameters used for this analysis.
170 char onePar[255] ;
171
172 sprintf(onePar,"--- AliAnaPhoton ---\n") ;
173 parList+=onePar ;
174 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
175 parList+=onePar ;
176 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
177 parList+=onePar ;
178 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
179 parList+=onePar ;
180 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
181 parList+=onePar ;
182 sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
183 parList+=onePar ;
184
185 //Get parameters set in base class.
186 parList += GetBaseParametersList() ;
187
188 //Get parameters set in PID class.
189 parList += GetCaloPID()->GetPIDParametersList() ;
190
191 //Get parameters set in FiducialCut class (not available yet)
192 //parlist += GetFidCut()->GetFidCutParametersList()
193
194 return new TObjString(parList) ;
195}
196
1c5acb87 197
198//________________________________________________________________________
199TList * AliAnaPhoton::GetCreateOutputObjects()
200{
477d6cee 201 // Create histograms to be saved in output file and
202 // store them in outputContainer
203 TList * outputContainer = new TList() ;
204 outputContainer->SetName("PhotonHistos") ;
4a745797 205
5a2dbc3c 206 Int_t nptbins = GetHistoPtBins();
207 Int_t nphibins = GetHistoPhiBins();
208 Int_t netabins = GetHistoEtaBins();
477d6cee 209 Float_t ptmax = GetHistoPtMax();
210 Float_t phimax = GetHistoPhiMax();
211 Float_t etamax = GetHistoEtaMax();
212 Float_t ptmin = GetHistoPtMin();
213 Float_t phimin = GetHistoPhiMin();
214 Float_t etamin = GetHistoEtaMin();
215
216 //Histograms of highest Photon identified in Event
217 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
218 fhPtPhoton->SetYTitle("N");
219 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
220 outputContainer->Add(fhPtPhoton) ;
221
222 fhPhiPhoton = new TH2F
223 ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
224 fhPhiPhoton->SetYTitle("#phi");
225 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
226 outputContainer->Add(fhPhiPhoton) ;
227
228 fhEtaPhoton = new TH2F
229 ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
230 fhEtaPhoton->SetYTitle("#eta");
231 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
232 outputContainer->Add(fhEtaPhoton) ;
233
234 if(IsDataMC()){
123fc3bd 235 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
236 fhDeltaE->SetXTitle("#Delta E (GeV)");
237 outputContainer->Add(fhDeltaE);
238
239 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
240 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
241 outputContainer->Add(fhDeltaPt);
242
243 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
244 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
245 outputContainer->Add(fhRatioE);
477d6cee 246
123fc3bd 247 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
248 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
249 outputContainer->Add(fhRatioPt);
250
251 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 252 fh2E->SetXTitle("E_{rec} (GeV)");
253 fh2E->SetYTitle("E_{gen} (GeV)");
123fc3bd 254 outputContainer->Add(fh2E);
255
256 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 257 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
258 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
123fc3bd 259 outputContainer->Add(fh2Pt);
260
591cc579 261 fhPtMCPhoton = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
262 fhPtMCPhoton->SetYTitle("N");
263 fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
264 outputContainer->Add(fhPtMCPhoton) ;
265
266 fhPhiMCPhoton = new TH2F
267 ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
268 fhPhiMCPhoton->SetYTitle("#phi");
269 fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
270 outputContainer->Add(fhPhiMCPhoton) ;
271
272 fhEtaMCPhoton = new TH2F
273 ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
274 fhEtaMCPhoton->SetYTitle("#eta");
275 fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
276 outputContainer->Add(fhEtaMCPhoton) ;
277
278 fhPtPrompt = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 279 fhPtPrompt->SetYTitle("N");
280 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
281 outputContainer->Add(fhPtPrompt) ;
282
283 fhPhiPrompt = new TH2F
a3aebfff 284 ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 285 fhPhiPrompt->SetYTitle("#phi");
286 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
287 outputContainer->Add(fhPhiPrompt) ;
288
289 fhEtaPrompt = new TH2F
a3aebfff 290 ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 291 fhEtaPrompt->SetYTitle("#eta");
292 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
293 outputContainer->Add(fhEtaPrompt) ;
294
591cc579 295 fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 296 fhPtFragmentation->SetYTitle("N");
297 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
298 outputContainer->Add(fhPtFragmentation) ;
299
300 fhPhiFragmentation = new TH2F
a3aebfff 301 ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 302 fhPhiFragmentation->SetYTitle("#phi");
303 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
304 outputContainer->Add(fhPhiFragmentation) ;
305
306 fhEtaFragmentation = new TH2F
a3aebfff 307 ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 308 fhEtaFragmentation->SetYTitle("#eta");
309 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
310 outputContainer->Add(fhEtaFragmentation) ;
311
a3aebfff 312 fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 313 fhPtISR->SetYTitle("N");
314 fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
315 outputContainer->Add(fhPtISR) ;
316
317 fhPhiISR = new TH2F
a3aebfff 318 ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 319 fhPhiISR->SetYTitle("#phi");
320 fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
321 outputContainer->Add(fhPhiISR) ;
322
323 fhEtaISR = new TH2F
a3aebfff 324 ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 325 fhEtaISR->SetYTitle("#eta");
326 fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
327 outputContainer->Add(fhEtaISR) ;
328
591cc579 329 fhPtPi0Decay = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 330 fhPtPi0Decay->SetYTitle("N");
331 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
332 outputContainer->Add(fhPtPi0Decay) ;
333
334 fhPhiPi0Decay = new TH2F
a3aebfff 335 ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 336 fhPhiPi0Decay->SetYTitle("#phi");
337 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
338 outputContainer->Add(fhPhiPi0Decay) ;
339
340 fhEtaPi0Decay = new TH2F
a3aebfff 341 ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 342 fhEtaPi0Decay->SetYTitle("#eta");
343 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
344 outputContainer->Add(fhEtaPi0Decay) ;
345
a3aebfff 346 fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 347 fhPtOtherDecay->SetYTitle("N");
348 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
349 outputContainer->Add(fhPtOtherDecay) ;
350
351 fhPhiOtherDecay = new TH2F
a3aebfff 352 ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 353 fhPhiOtherDecay->SetYTitle("#phi");
354 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
355 outputContainer->Add(fhPhiOtherDecay) ;
356
357 fhEtaOtherDecay = new TH2F
a3aebfff 358 ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 359 fhEtaOtherDecay->SetYTitle("#eta");
360 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
361 outputContainer->Add(fhEtaOtherDecay) ;
362
a3aebfff 363 fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 364 fhPtConversion->SetYTitle("N");
365 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
366 outputContainer->Add(fhPtConversion) ;
367
368 fhPhiConversion = new TH2F
a3aebfff 369 ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 370 fhPhiConversion->SetYTitle("#phi");
371 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
372 outputContainer->Add(fhPhiConversion) ;
373
374 fhEtaConversion = new TH2F
a3aebfff 375 ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 376 fhEtaConversion->SetYTitle("#eta");
377 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
378 outputContainer->Add(fhEtaConversion) ;
379
a3aebfff 380 fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 381 fhPtUnknown->SetYTitle("N");
382 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
383 outputContainer->Add(fhPtUnknown) ;
384
385 fhPhiUnknown = new TH2F
a3aebfff 386 ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 387 fhPhiUnknown->SetYTitle("#phi");
388 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
389 outputContainer->Add(fhPhiUnknown) ;
390
391 fhEtaUnknown = new TH2F
a3aebfff 392 ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 393 fhEtaUnknown->SetYTitle("#eta");
394 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
395 outputContainer->Add(fhEtaUnknown) ;
a3aebfff 396
477d6cee 397 }//Histos with MC
0c1383b5 398
477d6cee 399 return outputContainer ;
400
1c5acb87 401}
402
6639984f 403//____________________________________________________________________________
404void AliAnaPhoton::Init()
405{
406
407 //Init
408 //Do some checks
1e86c71e 409 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 410 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 411 abort();
412 }
1e86c71e 413 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 414 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 415 abort();
416 }
417
418}
419
420
1c5acb87 421//____________________________________________________________________________
422void AliAnaPhoton::InitParameters()
423{
424
425 //Initialize the parameters of the analysis.
a3aebfff 426 AddToHistogramsName("AnaPhoton_");
427
1c5acb87 428 fCalorimeter = "PHOS" ;
429 fMinDist = 2.;
430 fMinDist2 = 4.;
431 fMinDist3 = 5.;
1e86c71e 432 fMassCut = 0.03; //30 MeV
433
4cf55759 434 fTimeCutMin = -1;
435 fTimeCutMax = 9999999;
2ac125bf 436 fNCellsCut = 0;
437
1e86c71e 438 fRejectTrackMatch = kTRUE ;
439 fCheckConversion = kFALSE;
440 fAddConvertedPairsToAOD = kFALSE;
441
1c5acb87 442}
443
444//__________________________________________________________________
445void AliAnaPhoton::MakeAnalysisFillAOD()
446{
477d6cee 447 //Do analysis and fill aods
448 //Search for photons in fCalorimeter
4a745797 449 TObjArray * pl = 0x0;
477d6cee 450
451 //Get vertex for photon momentum calculation
233e0df8 452 Double_t vertex[] = {0,0,0} ; //vertex
453 Double_t vertex2[] = {0,0,0} ; //vertex from second input aod
edd59991 454 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
233e0df8 455 {
456 GetReader()->GetVertex(vertex);
457 if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
458 }
edd59991 459 //printf("Vertex 0: %f,%f,%f\n",vertex[0],vertex[1],vertex[2]);
460 //printf("Vertex 1: %f,%f,%f\n",vertex2[0],vertex2[1],vertex2[2]);
477d6cee 461 //Select the Calorimeter of the photon
462 if(fCalorimeter == "PHOS")
463 pl = GetAODPHOS();
464 else if (fCalorimeter == "EMCAL")
465 pl = GetAODEMCAL();
466
467 //Fill AODCaloClusters and AODParticle with PHOS aods
1e86c71e 468 TLorentzVector mom, mom2 ;
469 Int_t nCaloClusters = pl->GetEntriesFast();
470 Bool_t * indexConverted = new Bool_t[nCaloClusters];
471 for (Int_t i = 0; i < nCaloClusters; i++) indexConverted[i] = kFALSE;
472
473 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
474
475 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
233e0df8 476
ff45398a 477 //Cluster selection, not charged, with photon id and in fiducial cut
233e0df8 478
479 //Input from second AOD?
480 Int_t input = 0;
481 if (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo) input = 1 ;
482 else if(fCalorimeter == "PHOS" && GetReader()->GetAODPHOSNormalInputEntries() <= icalo) input = 1;
483
477d6cee 484 //Get Momentum vector,
1cd71065 485 if (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
486 else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
487
477d6cee 488 //If too small or big pt, skip it
489 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
4cf55759 490 Double_t tof = calo->GetTOF()*1e9;
2ac125bf 491
4cf55759 492 if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
2ac125bf 493
494 if(calo->GetNCells() <= fNCellsCut) continue;
495
591cc579 496 //printf("AliAnaPhoton::Current Event %d; Current File Name : %s, E %f, pT %f, Ecl %f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), mom.E(), mom.Pt(),calo->E());
497
477d6cee 498 //Check acceptance selection
ff45398a 499 if(IsFiducialCutOn()){
500 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 501 if(! in ) continue ;
502 }
591cc579 503
591cc579 504 //Create AOD for analysis
477d6cee 505 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
591cc579 506 Int_t label = calo->GetLabel(0);
591cc579 507 aodph.SetLabel(label);
233e0df8 508 aodph.SetInputFileIndex(input);
7cd4e982 509
477d6cee 510 //printf("Index %d, Id %d\n",icalo, calo->GetID());
511 //Set the indeces of the original caloclusters
512 aodph.SetCaloLabel(calo->GetID(),-1);
513 aodph.SetDetector(fCalorimeter);
514 if(GetDebug() > 1)
ff45398a 515 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
477d6cee 516
517 //Check Distance to Bad channel, set bit.
518 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
519 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
520 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
521 continue ;
522
a3aebfff 523 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
477d6cee 524
525 if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
526 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
527 else aodph.SetDistToBad(0) ;
528
ece3ffc2 529 //Skip matched clusters with tracks
530 if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
531
477d6cee 532 //Check PID
533 //PID selection or bit setting
534 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
535 //Get most probable PID, check PID weights (in MC this option is mandatory)
536 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
a3aebfff 537 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 538 //If primary is not photon, skip it.
539 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
540 }
541 else if(IsCaloPIDOn()){
ece3ffc2 542
477d6cee 543 //Get most probable PID, 2 options check PID weights
544 //or redo PID, recommended option for EMCal.
545 if(!IsCaloPIDRecalculationOn())
546 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
547 else
548 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
549
a3aebfff 550 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 551
552 //If cluster does not pass pid, not photon, skip it.
553 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
554
555 }
556 else{
557 //Set PID bits for later selection (AliAnaPi0 for example)
558 //GetPDG already called in SetPIDBits.
559 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
a3aebfff 560 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 561 }
562
a3aebfff 563 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
477d6cee 564
565 //Play with the MC stack if available
566 //Check origin of the candidates
567 if(IsDataMC()){
591cc579 568
78219bac 569 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabel(),GetReader(), aodph.GetInputFileIndex()));
591cc579 570 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 571 }//Work with stack also
572
1e86c71e 573
574 // Check if cluster comes from a conversion in the material in front of the calorimeter
575 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
576
577 if(fCheckConversion && nCaloClusters > 1){
578 Bool_t bConverted = kFALSE;
579 Int_t id2 = -1;
580
581 //Check if set previously as converted couple, if so skip its use.
582 if (indexConverted[icalo]) continue;
583
584 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
585 //Check if set previously as converted couple, if so skip its use.
586 if (indexConverted[jcalo]) continue;
587 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
588 AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics
589 calo2->GetMomentum(mom2,vertex);
590 //Check only certain regions
591 Bool_t in2 = kTRUE;
592 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
593 if(!in2) continue;
594
595 //Get mass of pair, if small, take this pair.
596 //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
597 if((mom+mom2).M() < fMassCut){
598 bConverted = kTRUE;
599 id2 = calo2->GetID();
600 indexConverted[jcalo]=kTRUE;
601 break;
602 }
603
604 }//Mass loop
605
606 if(bConverted){
607 if(fAddConvertedPairsToAOD){
608 //Create AOD of pair analysis
609 TLorentzVector mpair = mom+mom2;
610 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
611 aodpair.SetLabel(aodph.GetLabel());
612 aodpair.SetInputFileIndex(input);
613
614 //printf("Index %d, Id %d\n",icalo, calo->GetID());
615 //Set the indeces of the original caloclusters
616 aodpair.SetCaloLabel(calo->GetID(),id2);
617 aodpair.SetDetector(fCalorimeter);
618 aodpair.SetPdg(aodph.GetPdg());
619 aodpair.SetTag(aodph.GetTag());
620
621 //Add AOD with pair object to aod branch
622 AddAODParticle(aodpair);
623 //printf("\t \t both added pair\n");
624 }
625
626 //Do not add the current calocluster
627 continue;
628 }//converted pair
629 }//check conversion
630 //printf("\t \t added single cluster %d\n",icalo);
631
477d6cee 632 //Add AOD with photon object to aod branch
633 AddAODParticle(aodph);
634
635 }//loop
636
4a745797 637 delete [] indexConverted;
638
a3aebfff 639 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs \n");
477d6cee 640
1c5acb87 641}
642
643//__________________________________________________________________
644void AliAnaPhoton::MakeAnalysisFillHistograms()
645{
123fc3bd 646 //Do analysis and fill histograms
647
591cc579 648 // Access MC information in stack if requested, check that it exists.
649 AliStack * stack = 0x0;
650 TParticle * primary = 0x0;
651 TClonesArray * mcparticles0 = 0x0;
652 TClonesArray * mcparticles1 = 0x0;
653 AliAODMCParticle * aodprimary = 0x0;
654 if(IsDataMC()){
655
656 if(GetReader()->ReadStack()){
657 stack = GetMCStack() ;
658 if(!stack) {
659 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
660 abort();
661 }
662
663 }
664 else if(GetReader()->ReadAODMCParticles()){
665
666 //Get the list of MC particles
667 mcparticles0 = GetReader()->GetAODMCParticles(0);
668 if(!mcparticles0 && GetDebug() > 0) {
669 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
670 }
671 if(GetReader()->GetSecondInputAODTree()){
672 mcparticles1 = GetReader()->GetAODMCParticles(1);
673 if(!mcparticles1 && GetDebug() > 0) {
674 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
675 }
676 }
677
678 }
679 }// is data and MC
680
123fc3bd 681 //Loop on stored AOD photons
682 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
683 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
684
685 for(Int_t iaod = 0; iaod < naod ; iaod++){
477d6cee 686 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
687 Int_t pdg = ph->GetPdg();
688
689 if(GetDebug() > 3)
a3aebfff 690 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
477d6cee 691
692 //If PID used, fill histos with photons in Calorimeter fCalorimeter
123fc3bd 693 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
477d6cee 694 if(ph->GetDetector() != fCalorimeter) continue;
695
591cc579 696 if(GetDebug() > 2)
a3aebfff 697 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
477d6cee 698
699 //Fill photon histograms
123fc3bd 700 Float_t ptcluster = ph->Pt();
477d6cee 701 Float_t phicluster = ph->Phi();
702 Float_t etacluster = ph->Eta();
123fc3bd 703 Float_t ecluster = ph->E();
704
705 fhPtPhoton ->Fill(ptcluster);
477d6cee 706 fhPhiPhoton ->Fill(ptcluster,phicluster);
707 fhEtaPhoton ->Fill(ptcluster,etacluster);
123fc3bd 708
591cc579 709 //Play with the MC data if available
477d6cee 710 if(IsDataMC()){
123fc3bd 711
477d6cee 712 Int_t tag =ph->GetTag();
591cc579 713
714 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
715 {
716 fhPtMCPhoton ->Fill(ptcluster);
717 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
718 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
719
720 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
721 {
722 fhPtConversion ->Fill(ptcluster);
723 fhPhiConversion ->Fill(ptcluster,phicluster);
724 fhEtaConversion ->Fill(ptcluster,etacluster);
725 }
726
727 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
728 fhPtPrompt ->Fill(ptcluster);
729 fhPhiPrompt ->Fill(ptcluster,phicluster);
730 fhEtaPrompt ->Fill(ptcluster,etacluster);
731 }
732 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
733 {
734 fhPtFragmentation ->Fill(ptcluster);
735 fhPhiFragmentation ->Fill(ptcluster,phicluster);
736 fhEtaFragmentation ->Fill(ptcluster,etacluster);
737 }
738 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
739 {
740 fhPtISR ->Fill(ptcluster);
741 fhPhiISR ->Fill(ptcluster,phicluster);
742 fhEtaISR ->Fill(ptcluster,etacluster);
743 }
744 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
745 {
746 fhPtPi0Decay ->Fill(ptcluster);
747 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
748 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
749 }
750 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
751 {
752 fhPtOtherDecay ->Fill(ptcluster);
753 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
754 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
755 }
756 }
477d6cee 757 else{
591cc579 758 fhPtUnknown ->Fill(ptcluster);
759 fhPhiUnknown ->Fill(ptcluster,phicluster);
760 fhEtaUnknown ->Fill(ptcluster,etacluster);
761
762// printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
763// ph->GetLabel(),ph->Pt());
764// for(Int_t i = 0; i < 20; i++) {
765// if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
766// }
767// printf("\n");
768
477d6cee 769 }
591cc579 770
123fc3bd 771
591cc579 772 // Access MC information in stack if requested, check that it exists.
123fc3bd 773 Int_t label =ph->GetLabel();
774 if(label < 0) {
8d4bd2e1 775 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
123fc3bd 776 continue;
777 }
778
591cc579 779 Float_t eprim = 0;
780 Float_t ptprim = 0;
781 if(GetReader()->ReadStack()){
782
783 if(label >= stack->GetNtrack()) {
8d4bd2e1 784 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
591cc579 785 continue ;
786 }
787
788 primary = stack->Particle(label);
789 if(!primary){
790 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
791 continue;
792 }
793 eprim = primary->Energy();
794 ptprim = primary->Pt();
123fc3bd 795
123fc3bd 796 }
591cc579 797 else if(GetReader()->ReadAODMCParticles()){
798 //Check which is the input
7cd4e982 799 if(ph->GetInputFileIndex() == 0){
591cc579 800 if(!mcparticles0) continue;
801 if(label >= mcparticles0->GetEntriesFast()) {
802 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
803 label, mcparticles0->GetEntriesFast());
804 continue ;
805 }
806 //Get the particle
807 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
808
809 }
810 else {//Second input
811 if(!mcparticles1) continue;
591cc579 812 if(label >= mcparticles1->GetEntriesFast()) {
813 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
814 label, mcparticles1->GetEntriesFast());
815 continue ;
816 }
817 //Get the particle
818 aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
819
820 }//second input
821
822 if(!aodprimary){
823 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
824 continue;
825 }
826
827 eprim = aodprimary->E();
828 ptprim = aodprimary->Pt();
123fc3bd 829
591cc579 830 }
831
41e886c8 832 fh2E ->Fill(ecluster, eprim);
833 fh2Pt ->Fill(ptcluster, ptprim);
591cc579 834 fhDeltaE ->Fill(eprim-ecluster);
835 fhDeltaPt->Fill(ptprim-ptcluster);
836 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
837 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
123fc3bd 838
477d6cee 839 }//Histograms with MC
840
841 }// aod loop
842
1c5acb87 843}
844
845
846//__________________________________________________________________
847void AliAnaPhoton::Print(const Option_t * opt) const
848{
477d6cee 849 //Print some relevant parameters set for the analysis
850
851 if(! opt)
852 return;
853
854 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
855 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 856
477d6cee 857 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
858 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
859 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
860 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 861 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 862 printf("Check Pair Conversion = %d\n",fCheckConversion);
863 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
864 printf("Conversion pair mass cut = %f\n",fMassCut);
4cf55759 865 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 866 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 867 printf(" \n") ;
1c5acb87 868
869}