]> 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
Added ptHard to the information in the AliAODMCHeader, filled in the AODHandler
[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>
30#include <TObjString.h>
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"
477d6cee 40#include "AliFidutialCut.h"
41#include "AliAODCaloCluster.h"
1c5acb87 42
43ClassImp(AliAnaPhoton)
44
45//____________________________________________________________________________
46 AliAnaPhoton::AliAnaPhoton() :
47 AliAnaPartCorrBaseClass(), fCalorimeter(""),
a3aebfff 48 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
49 fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
7175a03a 50 //MC
123fc3bd 51 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
1c5acb87 52 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
53 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
6639984f 54 fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
1c5acb87 55 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
56 fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
57 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
58 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
59{
60 //default ctor
61
62 //Initialize parameters
63 InitParameters();
64
65}
66
67//____________________________________________________________________________
68AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
69 AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
70 fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
a3aebfff 71 fRejectTrackMatch(g.fRejectTrackMatch),
1c5acb87 72 fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
73 //MC
123fc3bd 74 fhDeltaE(g.fhDeltaE), fhDeltaPt(g.fhDeltaPt),
75 fhRatioE(g.fhRatioE), fhRatioPt(g.fhRatioPt),
76 fh2E(g.fh2E), fh2Pt(g.fh2Pt),
1c5acb87 77 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
78 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
6639984f 79 fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
1c5acb87 80 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
81 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
82 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
6639984f 83 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
84
1c5acb87 85{
86 // cpy ctor
87
88}
89
90//_________________________________________________________________________
91AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
92{
93 // assignment operator
94
95 if(&g == this) return *this;
96
97 fCalorimeter = g.fCalorimeter ;
98 fMinDist = g.fMinDist;
99 fMinDist2 = g.fMinDist2;
100 fMinDist3 = g.fMinDist3;
a3aebfff 101 fRejectTrackMatch = g.fRejectTrackMatch;
102
6639984f 103 fhPtPhoton = g.fhPtPhoton ;
1c5acb87 104 fhPhiPhoton = g.fhPhiPhoton ;
105 fhEtaPhoton = g.fhEtaPhoton ;
123fc3bd 106
107 fhDeltaE = g.fhDeltaE;
108 fhDeltaPt = g.fhDeltaPt;
109 fhRatioE = g.fhRatioE;
110 fhRatioPt = g.fhRatioPt;
111 fh2E = g.fh2E;
112 fh2Pt = g.fh2Pt;
1c5acb87 113
6639984f 114 fhPtPrompt = g.fhPtPrompt;
1c5acb87 115 fhPhiPrompt = g.fhPhiPrompt;
116 fhEtaPrompt = g.fhEtaPrompt;
6639984f 117 fhPtFragmentation = g.fhPtFragmentation;
1c5acb87 118 fhPhiFragmentation = g.fhPhiFragmentation;
119 fhEtaFragmentation = g.fhEtaFragmentation;
6639984f 120 fhPtISR = g.fhPtISR;
121 fhPhiISR = g.fhPhiISR;
122 fhEtaISR = g.fhEtaISR;
123 fhPtPi0Decay = g.fhPtPi0Decay;
1c5acb87 124 fhPhiPi0Decay = g.fhPhiPi0Decay;
125 fhEtaPi0Decay = g.fhEtaPi0Decay;
6639984f 126 fhPtOtherDecay = g.fhPtOtherDecay;
1c5acb87 127 fhPhiOtherDecay = g.fhPhiOtherDecay;
128 fhEtaOtherDecay = g.fhEtaOtherDecay;
6639984f 129 fhPtConversion = g. fhPtConversion;
1c5acb87 130 fhPhiConversion = g.fhPhiConversion;
131 fhEtaConversion = g.fhEtaConversion;
6639984f 132 fhPtUnknown = g.fhPtUnknown;
1c5acb87 133 fhPhiUnknown = g.fhPhiUnknown;
134 fhEtaUnknown = g.fhEtaUnknown;
135
136 return *this;
137
138}
139
140//____________________________________________________________________________
141AliAnaPhoton::~AliAnaPhoton()
142{
143 //dtor
144
145}
146
147
148//________________________________________________________________________
149TList * AliAnaPhoton::GetCreateOutputObjects()
150{
477d6cee 151 // Create histograms to be saved in output file and
152 // store them in outputContainer
153 TList * outputContainer = new TList() ;
154 outputContainer->SetName("PhotonHistos") ;
155
156 Int_t nptbins = GetHistoNPtBins();
157 Int_t nphibins = GetHistoNPhiBins();
158 Int_t netabins = GetHistoNEtaBins();
159 Float_t ptmax = GetHistoPtMax();
160 Float_t phimax = GetHistoPhiMax();
161 Float_t etamax = GetHistoEtaMax();
162 Float_t ptmin = GetHistoPtMin();
163 Float_t phimin = GetHistoPhiMin();
164 Float_t etamin = GetHistoEtaMin();
165
166 //Histograms of highest Photon identified in Event
167 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
168 fhPtPhoton->SetYTitle("N");
169 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
170 outputContainer->Add(fhPtPhoton) ;
171
172 fhPhiPhoton = new TH2F
173 ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
174 fhPhiPhoton->SetYTitle("#phi");
175 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
176 outputContainer->Add(fhPhiPhoton) ;
177
178 fhEtaPhoton = new TH2F
179 ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
180 fhEtaPhoton->SetYTitle("#eta");
181 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
182 outputContainer->Add(fhEtaPhoton) ;
183
184 if(IsDataMC()){
123fc3bd 185 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
186 fhDeltaE->SetXTitle("#Delta E (GeV)");
187 outputContainer->Add(fhDeltaE);
188
189 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
190 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
191 outputContainer->Add(fhDeltaPt);
192
193 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
194 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
195 outputContainer->Add(fhRatioE);
477d6cee 196
123fc3bd 197 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
198 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
199 outputContainer->Add(fhRatioPt);
200
201 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
202 fh2E->SetYTitle("E_{rec} (GeV)");
203 fh2E->SetXTitle("E_{gen} (GeV)");
204 outputContainer->Add(fh2E);
205
206 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
207 fh2Pt->SetYTitle("p_{T,rec} (GeV/c)");
208 fh2Pt->SetXTitle("p_{T,gen} (GeV/c)");
209 outputContainer->Add(fh2Pt);
210
a3aebfff 211 fhPtPrompt = new TH1F("hPtMCPrompt","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 212 fhPtPrompt->SetYTitle("N");
213 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
214 outputContainer->Add(fhPtPrompt) ;
215
216 fhPhiPrompt = new TH2F
a3aebfff 217 ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 218 fhPhiPrompt->SetYTitle("#phi");
219 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
220 outputContainer->Add(fhPhiPrompt) ;
221
222 fhEtaPrompt = new TH2F
a3aebfff 223 ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 224 fhEtaPrompt->SetYTitle("#eta");
225 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
226 outputContainer->Add(fhEtaPrompt) ;
227
a3aebfff 228 fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 229 fhPtFragmentation->SetYTitle("N");
230 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
231 outputContainer->Add(fhPtFragmentation) ;
232
233 fhPhiFragmentation = new TH2F
a3aebfff 234 ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 235 fhPhiFragmentation->SetYTitle("#phi");
236 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
237 outputContainer->Add(fhPhiFragmentation) ;
238
239 fhEtaFragmentation = new TH2F
a3aebfff 240 ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 241 fhEtaFragmentation->SetYTitle("#eta");
242 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
243 outputContainer->Add(fhEtaFragmentation) ;
244
a3aebfff 245 fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 246 fhPtISR->SetYTitle("N");
247 fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
248 outputContainer->Add(fhPtISR) ;
249
250 fhPhiISR = new TH2F
a3aebfff 251 ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 252 fhPhiISR->SetYTitle("#phi");
253 fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
254 outputContainer->Add(fhPhiISR) ;
255
256 fhEtaISR = new TH2F
a3aebfff 257 ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 258 fhEtaISR->SetYTitle("#eta");
259 fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
260 outputContainer->Add(fhEtaISR) ;
261
262 fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
263 fhPtPi0Decay->SetYTitle("N");
264 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
265 outputContainer->Add(fhPtPi0Decay) ;
266
267 fhPhiPi0Decay = new TH2F
a3aebfff 268 ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 269 fhPhiPi0Decay->SetYTitle("#phi");
270 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
271 outputContainer->Add(fhPhiPi0Decay) ;
272
273 fhEtaPi0Decay = new TH2F
a3aebfff 274 ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 275 fhEtaPi0Decay->SetYTitle("#eta");
276 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
277 outputContainer->Add(fhEtaPi0Decay) ;
278
a3aebfff 279 fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 280 fhPtOtherDecay->SetYTitle("N");
281 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
282 outputContainer->Add(fhPtOtherDecay) ;
283
284 fhPhiOtherDecay = new TH2F
a3aebfff 285 ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 286 fhPhiOtherDecay->SetYTitle("#phi");
287 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
288 outputContainer->Add(fhPhiOtherDecay) ;
289
290 fhEtaOtherDecay = new TH2F
a3aebfff 291 ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 292 fhEtaOtherDecay->SetYTitle("#eta");
293 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
294 outputContainer->Add(fhEtaOtherDecay) ;
295
a3aebfff 296 fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 297 fhPtConversion->SetYTitle("N");
298 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
299 outputContainer->Add(fhPtConversion) ;
300
301 fhPhiConversion = new TH2F
a3aebfff 302 ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 303 fhPhiConversion->SetYTitle("#phi");
304 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
305 outputContainer->Add(fhPhiConversion) ;
306
307 fhEtaConversion = new TH2F
a3aebfff 308 ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 309 fhEtaConversion->SetYTitle("#eta");
310 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
311 outputContainer->Add(fhEtaConversion) ;
312
a3aebfff 313 fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 314 fhPtUnknown->SetYTitle("N");
315 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
316 outputContainer->Add(fhPtUnknown) ;
317
318 fhPhiUnknown = new TH2F
a3aebfff 319 ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 320 fhPhiUnknown->SetYTitle("#phi");
321 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
322 outputContainer->Add(fhPhiUnknown) ;
323
324 fhEtaUnknown = new TH2F
a3aebfff 325 ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 326 fhEtaUnknown->SetYTitle("#eta");
327 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
328 outputContainer->Add(fhEtaUnknown) ;
a3aebfff 329
477d6cee 330 }//Histos with MC
331
332 //Save parameters used for analysis
333 TString parList ; //this will be list of parameters used for this analysis.
334 char onePar[255] ;
335
336 sprintf(onePar,"--- AliAnaPhoton ---\n") ;
337 parList+=onePar ;
338 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
339 parList+=onePar ;
340 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
341 parList+=onePar ;
342 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
343 parList+=onePar ;
344 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
345 parList+=onePar ;
a3aebfff 346 sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
347 parList+=onePar ;
477d6cee 348
349 //Get parameters set in base class.
350 parList += GetBaseParametersList() ;
351
352 //Get parameters set in PID class.
353 parList += GetCaloPID()->GetPIDParametersList() ;
354
355 //Get parameters set in FidutialCut class (not available yet)
356 //parlist += GetFidCut()->GetFidCutParametersList()
357
358 TObjString *oString= new TObjString(parList) ;
359 outputContainer->Add(oString);
360
361 return outputContainer ;
362
1c5acb87 363}
364
6639984f 365//____________________________________________________________________________
366void AliAnaPhoton::Init()
367{
368
369 //Init
370 //Do some checks
371 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
a3aebfff 372 printf("AliAnaPhoton::Init() - !!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 373 abort();
374 }
375 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
a3aebfff 376 printf("AliAnaPhoton::Init() - !!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 377 abort();
378 }
379
380}
381
382
1c5acb87 383//____________________________________________________________________________
384void AliAnaPhoton::InitParameters()
385{
386
387 //Initialize the parameters of the analysis.
388 SetOutputAODClassName("AliAODPWG4Particle");
a3aebfff 389 SetOutputAODName("PWG4Particle");
390
391 AddToHistogramsName("AnaPhoton_");
392
1c5acb87 393 fCalorimeter = "PHOS" ;
394 fMinDist = 2.;
395 fMinDist2 = 4.;
396 fMinDist3 = 5.;
a3aebfff 397 fRejectTrackMatch = kTRUE ;
1c5acb87 398}
399
400//__________________________________________________________________
401void AliAnaPhoton::MakeAnalysisFillAOD()
402{
477d6cee 403 //Do analysis and fill aods
404 //Search for photons in fCalorimeter
405 TRefArray * pl = new TRefArray();
406
407 //Get vertex for photon momentum calculation
408 Double_t vertex[]={0,0,0} ; //vertex ;
409 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
410 //Select the Calorimeter of the photon
411 if(fCalorimeter == "PHOS")
412 pl = GetAODPHOS();
413 else if (fCalorimeter == "EMCAL")
414 pl = GetAODEMCAL();
415
416 //Fill AODCaloClusters and AODParticle with PHOS aods
417 TLorentzVector mom ;
418
419 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
420 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
421 //Cluster selection, not charged, with photon id and in fidutial cut
422 //Get Momentum vector,
423 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
424 //If too small or big pt, skip it
425 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
426 //Check acceptance selection
427 if(IsFidutialCutOn()){
428 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
429 if(! in ) continue ;
430 }
431
432 //Create AOD for analysis
433 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
434 aodph.SetLabel(calo->GetLabel(0));
435 //printf("Index %d, Id %d\n",icalo, calo->GetID());
436 //Set the indeces of the original caloclusters
437 aodph.SetCaloLabel(calo->GetID(),-1);
438 aodph.SetDetector(fCalorimeter);
439 if(GetDebug() > 1)
a3aebfff 440 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
477d6cee 441
442 //Check Distance to Bad channel, set bit.
443 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
444 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
445 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
446 continue ;
447
a3aebfff 448 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
477d6cee 449
450 if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
451 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
452 else aodph.SetDistToBad(0) ;
453
454 //Check PID
455 //PID selection or bit setting
456 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
457 //Get most probable PID, check PID weights (in MC this option is mandatory)
458 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
a3aebfff 459 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 460 //If primary is not photon, skip it.
461 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
462 }
463 else if(IsCaloPIDOn()){
464 //Skip matched clusters with tracks
a3aebfff 465 if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
477d6cee 466
467 //Get most probable PID, 2 options check PID weights
468 //or redo PID, recommended option for EMCal.
469 if(!IsCaloPIDRecalculationOn())
470 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
471 else
472 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
473
a3aebfff 474 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 475
476 //If cluster does not pass pid, not photon, skip it.
477 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
478
479 }
480 else{
481 //Set PID bits for later selection (AliAnaPi0 for example)
482 //GetPDG already called in SetPIDBits.
483 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
a3aebfff 484 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 485 }
486
a3aebfff 487 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
477d6cee 488
489 //Play with the MC stack if available
490 //Check origin of the candidates
491 if(IsDataMC()){
492 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
a3aebfff 493 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate %d\n",aodph.GetTag());
477d6cee 494 }//Work with stack also
495
496 //Add AOD with photon object to aod branch
497 AddAODParticle(aodph);
498
499 }//loop
500
a3aebfff 501 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs \n");
477d6cee 502
1c5acb87 503}
504
505//__________________________________________________________________
506void AliAnaPhoton::MakeAnalysisFillHistograms()
507{
123fc3bd 508 //Do analysis and fill histograms
509
510 // Access MC information in stack if requested, check that it exists.
511 AliStack * stack = 0x0;
512 TParticle * primary = 0x0;
513 if(IsDataMC()){
514 stack = GetMCStack() ;
515 if(!stack) {
516 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
517 abort();
518 }
519 }
520
521 //Loop on stored AOD photons
522 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
523 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
524
525 for(Int_t iaod = 0; iaod < naod ; iaod++){
477d6cee 526 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
527 Int_t pdg = ph->GetPdg();
528
529 if(GetDebug() > 3)
a3aebfff 530 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
477d6cee 531
532 //If PID used, fill histos with photons in Calorimeter fCalorimeter
123fc3bd 533 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
477d6cee 534 if(ph->GetDetector() != fCalorimeter) continue;
535
536 if(GetDebug() > 1)
a3aebfff 537 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
477d6cee 538
539 //Fill photon histograms
123fc3bd 540 Float_t ptcluster = ph->Pt();
477d6cee 541 Float_t phicluster = ph->Phi();
542 Float_t etacluster = ph->Eta();
123fc3bd 543 Float_t ecluster = ph->E();
544
545 fhPtPhoton ->Fill(ptcluster);
477d6cee 546 fhPhiPhoton ->Fill(ptcluster,phicluster);
547 fhEtaPhoton ->Fill(ptcluster,etacluster);
123fc3bd 548
549 //Play with the MC stack if available
477d6cee 550 if(IsDataMC()){
123fc3bd 551
477d6cee 552 Int_t tag =ph->GetTag();
553 if(tag == AliMCAnalysisUtils::kMCPrompt){
123fc3bd 554 fhPtPrompt ->Fill(ptcluster);
555 fhPhiPrompt ->Fill(ptcluster,phicluster);
556 fhEtaPrompt ->Fill(ptcluster,etacluster);
477d6cee 557 }
558 else if(tag==AliMCAnalysisUtils::kMCFragmentation)
559 {
560 fhPtFragmentation ->Fill(ptcluster);
561 fhPhiFragmentation ->Fill(ptcluster,phicluster);
562 fhEtaFragmentation ->Fill(ptcluster,etacluster);
563 }
564 else if(tag==AliMCAnalysisUtils::kMCISR)
565 {
566 fhPtISR ->Fill(ptcluster);
567 fhPhiISR ->Fill(ptcluster,phicluster);
568 fhEtaISR ->Fill(ptcluster,etacluster);
569 }
570 else if(tag==AliMCAnalysisUtils::kMCPi0Decay)
571 {
572 fhPtPi0Decay ->Fill(ptcluster);
573 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
574 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
575 }
576 else if(tag==AliMCAnalysisUtils::kMCEtaDecay || tag==AliMCAnalysisUtils::kMCOtherDecay)
577 {
578 fhPtOtherDecay ->Fill(ptcluster);
579 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
580 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
581 }
582 else if(tag==AliMCAnalysisUtils::kMCConversion)
583 {
584 fhPtConversion ->Fill(ptcluster);
585 fhPhiConversion ->Fill(ptcluster,phicluster);
586 fhEtaConversion ->Fill(ptcluster,etacluster);
587
588 }
589 else{
590 fhPtUnknown ->Fill(ptcluster);
591 fhPhiUnknown ->Fill(ptcluster,phicluster);
592 fhEtaUnknown ->Fill(ptcluster,etacluster);
593 }
123fc3bd 594
595 Int_t label =ph->GetLabel();
596 if(label < 0) {
597 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
598 continue;
599 }
600
601 if(label >= stack->GetNtrack()) {
602 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
603 continue ;
604 }
605
606 primary = stack->Particle(label);
607 if(!primary){
608 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
609 continue;
610 }
611
612 fh2E->Fill(primary->Energy(),ecluster);
613 fh2Pt->Fill(primary->Pt(), ptcluster);
614 fhDeltaE->Fill(primary->Energy()-ecluster);
615 fhDeltaPt->Fill(primary->Pt()-ptcluster);
616 if(primary->Energy()!=0) fhRatioE ->Fill(ecluster/primary->Energy());
617 if(primary->Pt()!=0) fhRatioPt ->Fill(ptcluster/primary->Pt());
618
477d6cee 619 }//Histograms with MC
620
621 }// aod loop
622
1c5acb87 623}
624
625
626//__________________________________________________________________
627void AliAnaPhoton::Print(const Option_t * opt) const
628{
477d6cee 629 //Print some relevant parameters set for the analysis
630
631 if(! opt)
632 return;
633
634 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
635 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 636
477d6cee 637 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
638 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
639 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
640 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 641 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
642
477d6cee 643 printf(" \n") ;
1c5acb87 644
645}