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