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