]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
changed addtask PbPb
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_Tagging / AliAnalysisTaskTaggedPhotons.cxx
CommitLineData
213dd723 1/**************************************************************************
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 is 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
16/* $Id: */
17
18//_________________________________________________________________________
19// Analysis for Tagged Photons
20// Prepares all necessary histograms for calculation of
21// the yield of pi0 decay photon in calorimeter:
22// Marks photons which makes pi0 with some other and
23// fill invariant mass distributions for estimate background below pi0
24// peak so that estimate proportion of fake pairs.
25// Fills as well controll MC histograms with clasification of the photon origin
26// and check of the proportion of truly tagged photons.
27//
28//
29//*-- Dmitry Blau, Dmitri Peresunko
30//////////////////////////////////////////////////////////////////////////////
31
32#include <TH1.h>
33#include <TH2.h>
34#include <TH3.h>
35#include <THashList.h>
36#include <TFile.h>
37#include <TROOT.h>
38
39#include "AliAnalysisTaskTaggedPhotons.h"
40#include "AliAnalysisManager.h"
41#include "AliAODEvent.h"
42#include "AliAODEvent.h"
43#include "AliVCluster.h"
44#include "AliAODPWG4Particle.h"
45#include "AliAnalysisManager.h"
46#include "AliLog.h"
47#include "TGeoManager.h"
48#include "AliMCAnalysisUtils.h"
49#include "AliMCEventHandler.h"
50#include "AliMCEvent.h"
51#include "AliStack.h"
52#include "AliPHOSGeometry.h"
53#include "AliTriggerAnalysis.h"
54#include "AliEMCALGeometry.h"
55
56
57ClassImp(AliAnalysisTaskTaggedPhotons)
58
59//______________________________________________________________________________
60AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() :
61 AliAnalysisTaskSE(),
62 fPHOSgeom(0x0),
63 fOutputContainer(0x0),
64 fStack(0x0),
65 fTrackEvent(0x0),
66 fPHOSEvent(0x0),
67 fCurrentMixedList(0x0),
68 fTriggerAnalysis(0x0),
69 fZmax(0.),
70 fZmin(0.),
71 fPhimax(0.),
72 fPhimin(0.),
73 fCentrality(0.),
74 fCentBin(0)
75{
76 //Deafult constructor
77 //no memory allocations
78 for(Int_t i=0;i<10;i++)
79 for(Int_t j=0;j<2;j++)
80 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
bffa2e93 81 for(Int_t mod=0; mod<6; mod++)
82 fPHOSBadMap[mod]=0x0 ;
213dd723 83}
84//______________________________________________________________________________
85AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) :
86 AliAnalysisTaskSE(name),
87 fPHOSgeom(0x0),
88 fOutputContainer(0x0),
89 fStack(0x0),
90 fTrackEvent(0x0),
91 fPHOSEvent(0x0),
92 fCurrentMixedList(0x0),
93 fTriggerAnalysis(new AliTriggerAnalysis),
94 fZmax(-60.),
95 fZmin(60.),
96 fPhimax(250.),
97 fPhimin(320.),
98 fCentrality(0.),
99 fCentBin(0)
100{
101 // Constructor.
102
103 // Output slots #0 write into a TH1 container
104 DefineOutput(1,THashList::Class());
105 // Set bad channel map (empty so far)
106 for(Int_t i=0;i<1;i++)
107 for(Int_t j=0;j<5;j++)
108 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
109 // Create AOD track cut
110
111 for(Int_t mod=0; mod<6; mod++)
112 fPHOSBadMap[mod]=0x0 ;
113
114
115}
116
117//____________________________________________________________________________
118AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :
119 AliAnalysisTaskSE(ap.GetName()),
120 fPHOSgeom(0x0),
121 fOutputContainer(0x0),
122 fStack(0x0),
123 fTrackEvent(0x0),
124 fPHOSEvent(0x0),
125 fCurrentMixedList(0x0),
126 fTriggerAnalysis(new AliTriggerAnalysis),
127 fZmax(-60.),
128 fZmin(60.),
129 fPhimax(250.),
130 fPhimin(320.),
131 fCentrality(0.),
132 fCentBin(0)
133{
134 // cpy ctor
135 fZmax=ap.fZmax ;
136 fZmin=ap.fZmin ;
137 fPhimax=ap.fPhimax ;
138 fPhimin=ap.fPhimin ;
139 for(Int_t i=0;i<1;i++)
140 for(Int_t j=0;j<5;j++)
141 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
bffa2e93 142 for(Int_t mod=0; mod<6; mod++)
143 fPHOSBadMap[mod]=0x0 ;
144
213dd723 145}
146
147//_____________________________________________________________________________
148AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)
149{
150// assignment operator
151
152 this->~AliAnalysisTaskTaggedPhotons();
153 new(this) AliAnalysisTaskTaggedPhotons(ap);
154 return *this;
155}
156
157//______________________________________________________________________________
158AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()
159{
160 // dtor
161 if(fOutputContainer && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
162 fOutputContainer->Clear() ;
163 delete fOutputContainer ;
164 }
165 for(Int_t i=0;i<10;i++)
166 for(Int_t j=0;j<2;j++)
167 if(fPHOSEvents[i][j]){
168 delete fPHOSEvents[i][j] ;
169 fPHOSEvents[i][j]=0x0 ;
170 }
171}
172//________________________________________________________________________
173void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
174{
175
176
177 // Create the outputs containers
178 if(fOutputContainer != NULL){
179 delete fOutputContainer;
180 }
181 fOutputContainer = new THashList();
182 fOutputContainer->SetOwner(kTRUE);
183 fOutputContainer->SetName(GetName()) ;
184
185 //QA histograms
186 fOutputContainer->Add(new TH1F("hSelEvents","Event selection", 10,0.,10.)) ;
187
188 //vertex distribution
189 fOutputContainer->Add(new TH1F("hNvertexTracks","N of primary tracks from the primary vertex",150,0.,150.));
190 fOutputContainer->Add(new TH1F("hZvertex","Z vertex",200,-50.,+50.));
191 fOutputContainer->Add(new TH2F("hTrackMult","Charged track multiplicity",100,0.,100.,250,0.,500.));
192
193 //centrality
194 fOutputContainer->Add(new TH1F("hCentrality","Ccentrality",100,0.,100.));
195
196 fOutputContainer->Add(new TH2F("hTOF","cluster TOF",200,0.,20.,300,-3.e-6,6.e-6));
197
198
199 fOutputContainer->Add(new TH2F("hCluNXZM1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
200 fOutputContainer->Add(new TH2F("hCluNXZM2","Clu (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
201 fOutputContainer->Add(new TH2F("hCluNXZM3","Clu (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
202 fOutputContainer->Add(new TH2F("hCluEXZM1","Clu E(X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
203 fOutputContainer->Add(new TH2F("hCluEXZM2","Clu E(X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
204 fOutputContainer->Add(new TH2F("hCluEXZM3","Clu E(X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
205
206 fOutputContainer->Add(new TH2F("hCluArea2M1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
207 fOutputContainer->Add(new TH2F("hCluArea2M2","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
208 fOutputContainer->Add(new TH2F("hCluArea2M3","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
209 fOutputContainer->Add(new TH2F("hCluArea3M1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
210 fOutputContainer->Add(new TH2F("hCluArea3M2","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
211 fOutputContainer->Add(new TH2F("hCluArea3M3","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
212
213 char cPID[4][5] ;
214 snprintf(cPID[0],5,"All") ;
215 snprintf(cPID[1],5,"Disp");
216 snprintf(cPID[2],5,"CPV") ;
217 snprintf(cPID[3],5,"Both");
218
219
220 const Int_t nPt=400 ;
221 const Double_t ptMax=40. ;
222 const Int_t nM=400 ;
223 const Double_t mMax=1. ;
224
225 const Int_t nCenBin=5;
226 for(Int_t cen=0; cen<nCenBin; cen++){
227
228 for(Int_t iPID=0; iPID<4; iPID++){
229
230 //Inclusive spectra
231 fOutputContainer->Add(new TH1F(Form("hPhot_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
232 fOutputContainer->Add(new TH1F(Form("hPhot_Dist2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
233 fOutputContainer->Add(new TH1F(Form("hPhot_Dist3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
234
235 fOutputContainer->Add(new TH1F(Form("hPhot_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
236 fOutputContainer->Add(new TH1F(Form("hPhot_Area2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
237 fOutputContainer->Add(new TH1F(Form("hPhot_Area3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
238
239 fOutputContainer->Add(new TH1F(Form("hPhot_nStrictTagged_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
240
241 for(Int_t itag=0; itag<18; itag++){
242 fOutputContainer->Add(new TH1F(Form("hPhot_nTagged%d_Area1_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
243 fOutputContainer->Add(new TH1F(Form("hPhot_nTagged%d_Area2_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
244 fOutputContainer->Add(new TH1F(Form("hPhot_nTagged%d_Area3_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
245 }
246 for(Int_t kind=1; kind<33; kind*=2){
247 fOutputContainer->Add(new TH1F(Form("hPhot_Isolation%d_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
248 fOutputContainer->Add(new TH1F(Form("hPhot_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
249 fOutputContainer->Add(new TH1F(Form("hPhot_nStrictTagged_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
250 fOutputContainer->Add(new TH1F(Form("hPhot_nTagged_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
251 }
252 }
253
254
255 fOutputContainer->Add(new TH1F(Form("hTaggedMult_cent%d",cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
256
257 //Invariant mass distributions for fake corrections
258
259 for(Int_t iPID=0; iPID<4; iPID++){
260 fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
261 fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
262 fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
263 fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin1_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
264 fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin2_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
265 fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin3_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
266
267 fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
268 fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
269 fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
270 fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin1_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
271 fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin2_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
272 fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin3_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
273 }
274
275
276 fOutputContainer->Add(new TH2F(Form("QA_Cone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
277 fOutputContainer->Add(new TH2F(Form("QA_Cone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
278 fOutputContainer->Add(new TH2F(Form("QA_Cone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
279 }//centrality
280
281 //MC
282 char partName[4][10] ;
283 snprintf(partName[0],10,"gamma") ;
284 snprintf(partName[1],10,"pi0");
285 snprintf(partName[2],10,"eta") ;
286 snprintf(partName[3],10,"omega");
287
288 if(AliAnalysisManager::GetAnalysisManager()){
289 AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
290 if(mctruth){
b5f58f9c 291
292 fOutputContainer->Add(new TH1F("hMCConversionRadius","Clusters without label",600,0.,600.)) ;
293 fOutputContainer->Add(new TH2F("hMCRecPi0Vtx","Secondary pi0s",100,0.,10.,600,0.,600.)) ;
294 fOutputContainer->Add(new TH2F("hMCRecEtaVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
295 fOutputContainer->Add(new TH2F("hMCRecOmegaVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
296 fOutputContainer->Add(new TH2F("hMCRecEtaprVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
297 fOutputContainer->Add(new TH2F("hMCRecK0sVtx","Secondary K0s",100,0.,10.,600,0.,600.)) ;
298 fOutputContainer->Add(new TH2F("hMCRecK0lVtx","Secondary K0l",100,0.,10.,600,0.,600.)) ;
299 fOutputContainer->Add(new TH2F("hMCGammaPi0MisConvR","Converted photons",400,0.,40.,600,0.,600.)) ;
300
301 for(Int_t cen=0; cen<nCenBin; cen++){
213dd723 302 for(Int_t ipart=0; ipart<4; ipart++){
303 fOutputContainer->Add(new TH2F(Form("hMC_%s_vertex",partName[ipart]),"vertex",nPt,0.,ptMax,150,0.,600.)) ;
304 fOutputContainer->Add(new TH1F(Form("hMC_all_%s",partName[ipart]),"Spectum (full rapifity)",nPt,0.,ptMax)) ;
305 fOutputContainer->Add(new TH1F(Form("hMC_unitEta_%s",partName[ipart]),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
306 fOutputContainer->Add(new TH1F(Form("hMC_rap_%s",partName[ipart]),"Rapidity",100,-5.,5.)) ;
307 fOutputContainer->Add(new TH1F(Form("hMC_phi_%s",partName[ipart]),"Azimuthal angle",100,0.,TMath::TwoPi())) ;
308 }
309
310
b5f58f9c 311 fOutputContainer->Add(new TH2F(Form("LabelsNPrim_cent%d",cen),"vertex",nPt,0.,ptMax,20,0.,20.)) ;
312 fOutputContainer->Add(new TH1F(Form("LabelsGamma_cent%d",cen),"vertex",nPt,0.,ptMax)) ;
313 fOutputContainer->Add(new TH2F(Form("LabelsGammaE_cent%d",cen),"vertex",nPt,0.,ptMax,100,0.,2.)) ;
213dd723 314
b5f58f9c 315 fOutputContainer->Add(new TH1F(Form("hMCRecNoLabel_cent%d",cen),"Clusters without label",nPt,0.,ptMax)) ;
316 fOutputContainer->Add(new TH1F(Form("hMCConversionRadius_cent%d",cen),"Clusters without label",600,0.,600.)) ;
317 fOutputContainer->Add(new TH2F(Form("hMCRecPi0Vtx_cent%d",cen),"Secondary pi0s",100,0.,10.,600,0.,600.)) ;
318 fOutputContainer->Add(new TH2F(Form("hMCRecEtaVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
319 fOutputContainer->Add(new TH2F(Form("hMCRecOmegaVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
320 fOutputContainer->Add(new TH2F(Form("hMCRecEtaprVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
321 fOutputContainer->Add(new TH2F(Form("hMCRecK0sVtx_cent%d",cen),"Secondary K0s",100,0.,10.,600,0.,600.)) ;
322 fOutputContainer->Add(new TH2F(Form("hMCRecK0lVtx_cent%d",cen),"Secondary K0l",100,0.,10.,600,0.,600.)) ;
323 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisConvR_cent%d",cen),"Converted photons",400,0.,40.,600,0.,600.)) ;
213dd723 324
325
b5f58f9c 326 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0PrimMgg_cent%d",cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
327 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0RecMgg_cent%d",cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
213dd723 328
329 for(Int_t iPID=0; iPID<4; iPID++){
b5f58f9c 330 fOutputContainer->Add(new TH1F(Form("hMCRecNoPrim_%s_cent%d",cPID[iPID],cen),"Reconstructed particles withour primary",nPt,0.,ptMax)) ;
331 fOutputContainer->Add(new TH1F(Form("hMCRecGamma_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: Gamma",nPt,0.,ptMax)) ;
332 fOutputContainer->Add(new TH1F(Form("hMCRecPhotConv_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: e+-",nPt,0.,ptMax)) ;
333 fOutputContainer->Add(new TH1F(Form("hMCRecPi0_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: pi0",nPt,0.,ptMax)) ;
334 fOutputContainer->Add(new TH1F(Form("hMCRecEta_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: eta",nPt,0.,ptMax)) ;
335 fOutputContainer->Add(new TH1F(Form("hMCRecOmega_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: Gamma",nPt,0.,ptMax)) ;
336 fOutputContainer->Add(new TH1F(Form("hMCRecEtapr_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: eta prime",nPt,0.,ptMax)) ;
337 fOutputContainer->Add(new TH1F(Form("hMCRecPbar_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: bar(p)",nPt,0.,ptMax)) ;
338 fOutputContainer->Add(new TH1F(Form("hMCRecNbar_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: bar(n)",nPt,0.,ptMax)) ;
339 fOutputContainer->Add(new TH1F(Form("hMCRecPipm_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: pipm",nPt,0.,ptMax)) ;
340 fOutputContainer->Add(new TH1F(Form("hMCRecN_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: n",nPt,0.,ptMax)) ;
341 fOutputContainer->Add(new TH1F(Form("hMCRecP_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: p",nPt,0.,ptMax)) ;
342 fOutputContainer->Add(new TH1F(Form("hMCRecKpm_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K+-",nPt,0.,ptMax)) ;
343 fOutputContainer->Add(new TH1F(Form("hMCRecK0s_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0s",nPt,0.,ptMax)) ;
344 fOutputContainer->Add(new TH1F(Form("hMCRecK0l_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
345 fOutputContainer->Add(new TH1F(Form("hMCRecUnknownCh_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
346 fOutputContainer->Add(new TH1F(Form("hMCRecUnknownNeu_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
213dd723 347
348 //Decay photons
b5f58f9c 349 fOutputContainer->Add(new TH1F(Form("hMCRecGammaDir_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, no primary",nPt,0.,ptMax)) ;
350 fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
351 fOutputContainer->Add(new TH1F(Form("hMCRecGammaEta_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from eta",nPt,0.,ptMax)) ;
352 fOutputContainer->Add(new TH1F(Form("hMCRecGammaOmega_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from omega",nPt,0.,ptMax)) ;
353 fOutputContainer->Add(new TH1F(Form("hMCRecGammaOther_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from other",nPt,0.,ptMax)) ;
213dd723 354
355 //Pi0 decay photons
b5f58f9c 356 fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0Dalitz_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
357 fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0NoStack_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
358 fOutputContainer->Add(new TH1F(Form("hMCGammaPi02Gamma_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
359 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0Rec_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
360 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0RecSoft_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
361 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0RecCells_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
362
363 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0Tagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
364 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0FakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
365 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0TrueTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
366 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MultyTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
367 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeo_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
368 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA1_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
369 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA2_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
370 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA3_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
371 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisConv_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
372 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisEmin_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
373 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisOther_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
374 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisFakePrim_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
375 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisFoundPrim_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
376 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisNPhot_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
213dd723 377
378 //all clusters fake tagged
b5f58f9c 379 fOutputContainer->Add(new TH1F(Form("hMCAllFakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
213dd723 380
b5f58f9c 381 fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
382 fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
383
213dd723 384 }
b5f58f9c 385 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisPartner_cent%d",cen),"Spectrum of missed partners",nPt,0.,ptMax)) ;
386 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisPartnerEtaPhi_cent%d",cen),"Spectrum of missed partners",100,-0.2,0.2,100,4.5,5.6)) ;
387 }
213dd723 388 }
389 }//If MC handler exists...
390
391 for(Int_t i=0;i<1;i++)
392 for(Int_t j=0;j<5;j++)
393 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
394
395
396 PostData(1, fOutputContainer);
397
398
399}
400
401//______________________________________________________________________________
402void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
403{
404 //Select events
405 //Select photons
406 //Fill QA histograms
407 //Fill Tagging histogsms
408
409 // Event selection flags
410
411 FillHistogram("hSelEvents",1) ;
412
b5f58f9c 413 AliVEvent* event = (AliVEvent*)InputEvent();
213dd723 414 if(!event){
415 AliDebug(1,"No event") ;
416 PostData(1, fOutputContainer);
417 return;
418 }
419 FillHistogram("hSelEvents",2) ;
420
421 //read geometry if not read yet
422 if(fPHOSgeom==0)
423 InitGeometry() ;
424
425 //MC stack init
426 fStack=0x0 ;
427 if(AliAnalysisManager::GetAnalysisManager()){
428 AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
429 if(mctruth)
430 fStack = mctruth->MCEvent()->Stack();
431 }
432
433
434 // Checks if we have a primary vertex
435 // Get primary vertices form AOD
436
437/*
438 if(event->GetPrimaryVertexTracks()->GetNContributors()>0)
439 eventVtxExist = kTRUE;
440 else
441 if(event->GetPrimaryVertexSPD()->GetNContributors()>0)
442 eventVtxExist = kTRUE;
443*/
213dd723 444
445 Double_t vtx5[3];
b5f58f9c 446 vtx5[0] = event->GetPrimaryVertex()->GetX();
447 vtx5[1] = event->GetPrimaryVertex()->GetY();
448 vtx5[2] = event->GetPrimaryVertex()->GetZ();
449
450 FillHistogram("hNvertexTracks",event->GetPrimaryVertex()->GetNContributors());
451 FillHistogram("hZvertex" ,vtx5[2]);
452 if (TMath::Abs(vtx5[2]) > 10. ){
213dd723 453 PostData(1, fOutputContainer);
454 return ;
455 }
456
457 FillHistogram("hSelEvents",3) ;
458
459
b5f58f9c 460// if (event->IsPileupFromSPD()){
461// PostData(1, fOutputContainer);
462// return ;
463// }
213dd723 464
465 FillHistogram("hSelEvents",4) ;
466
467 //centrality
468 AliCentrality *centrality = event->GetCentrality();
469 if( centrality )
470 fCentrality=centrality->GetCentralityPercentile("V0M");
471 else {
472 AliError("Event has 0x0 centrality");
473 fCentrality = -1.;
474 }
475 FillHistogram("hCentrality",fCentrality) ;
476
477 if(fCentrality<0. || fCentrality>=100.){
478 PostData(1, fOutputContainer);
479 return ;
480 }
481 fCentBin = (Int_t)(fCentrality/20.) ;
482
483 FillHistogram("hSelEvents",5) ;
484
485
486 //Vtx class z-bin
487 Int_t zvtx = 0 ;
488
489 //Calculate charged multiplicity
490 Int_t trackMult = 0;
491 if(fTrackEvent)
492 fTrackEvent->Clear() ;
493 else
494 fTrackEvent = new TClonesArray("AliAODPWG4Particle",event->GetNumberOfTracks()) ;
495 for (Int_t i=0;i<event->GetNumberOfTracks();++i) {
b5f58f9c 496 AliVParticle *track = event->GetTrack(i) ;
213dd723 497 if(TMath::Abs(track->Eta())< 0.9){
498 if(trackMult>=fTrackEvent->GetSize())
499 fTrackEvent->Expand(2*trackMult) ;
500 new ((*fTrackEvent)[trackMult]) AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),track->P());
501 trackMult++;
502 }
213dd723 503 }
504 FillHistogram("hTrackMult",fCentrality,trackMult+0.5) ;
505
506 if(!fPHOSEvents[zvtx][fCentBin])
507 fPHOSEvents[zvtx][fCentBin]=new TList() ;
508 fCurrentMixedList = fPHOSEvents[zvtx][fCentBin] ;
509
510 const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
511
512 //---------Select photons-------------------
513 Int_t multClust = event->GetNumberOfCaloClusters();
514 if(!fPHOSEvent)
515 fPHOSEvent = new TClonesArray("AliAODPWG4Particle",multClust);
516 else
517 fPHOSEvent->Clear() ;
518 Int_t inList = 0; //counter of caloClusters
519
520 for (Int_t i=0; i<multClust; i++) {
b5f58f9c 521 AliVCluster * clu = event->GetCaloCluster(i);
213dd723 522
523 if(!clu->IsPHOS())
524 continue ;
525
526 if(clu->E()<0.1)
527 continue;
528
529 if(clu->GetNCells()<3)
530 continue ;
531
532 if(clu->GetM02()<0.2)
533 continue ;
534
535 FillHistogram("hTOF",clu->E(),clu->GetTOF()) ;
b5f58f9c 536 if(TMath::Abs(clu->GetTOF())>100.e-9)
213dd723 537 continue ;
538
539 if(clu->GetDistanceToBadChannel()<2.5)
540 continue ;
541
542 Float_t pos[3] ;
543 clu->GetPosition(pos) ;
544 Int_t fidArea=GetFiducialArea(pos) ;
545 if(fidArea==0) //Bad cell
546 continue;
547
548
549 TVector3 global1(pos) ;
550 Int_t relId[4] ;
551 fPHOSgeom->GlobalPos2RelId(global1,relId) ;
552 Int_t mod = relId[0] ;
553 Int_t cellX = relId[2];
554 Int_t cellZ = relId[3] ;
555
556 FillHistogram(Form("hCluNXZM%d",mod),cellX,cellZ,1.);
557 FillHistogram(Form("hCluEXZM%d",mod),cellX,cellZ,clu->E());
558 if(fidArea>1){
559 FillHistogram(Form("hCluArea2M%d",mod),cellX,cellZ,1.);
560 if(fidArea>2){
561 FillHistogram(Form("hCluArea3M%d",mod),cellX,cellZ,1.);
562 }
563 }
564
565 TLorentzVector momentum ;
566 clu->GetMomentum(momentum, vtx5);
567 AliAODPWG4Particle *p = new ((*fPHOSEvent)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),clu->E() );
568 inList++;
569
570 Int_t isolation = EvalIsolation(&momentum) ;
571 p->SetBtag(isolation) ;
572
573 p->SetDistToBad((Int_t)(1.+clu->GetDistanceToBadChannel()/2.2));
574
575 p->SetTag(0); //Strict PID pi0 partner not found
576 p->SetTagged(kFALSE); //Reconstructed pairs found
577 Bool_t sure = 0;
578 p->SetFiducialArea(fidArea) ;
579
580 if(fStack){
581 //This is primary entered PHOS
582 Int_t primLabel=FindPrimary(clu,sure) ;
583 //Look what particle left vertex
584 if(primLabel>-1){
585 TParticle * prim = fStack->Particle(primLabel) ;
586 Int_t iparent=primLabel;
587 TParticle * parent = prim;
588 while((parent->R() > rcut) && (iparent>-1)){
589 iparent=parent->GetFirstMother();
590 parent=fStack->Particle(iparent);
591 }
592 p->SetCaloLabel(primLabel,iparent); //This and partner cluster
593 }
594 else{
595 p->SetCaloLabel(-1,-1); //This and partner cluster
596 }
597
598 }
599 else{
600 //This is primary at vertex(R<1cm)
601 p->SetCaloLabel(-1,-1); //This and partner cluster
602 }
603 p->SetLabel(i); //Cluster index
604 //PID criteria
605// p->SetDispBit(clu->Chi2()) ;
606 p->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
607 p->SetTOFBit(TestTOF(clu->GetTOF(),clu->E())) ;
608 p->SetChargedBit(clu->GetEmcCpvDistance()>10.) ;
609 }
610
611 FillMCHistos() ;
612 FillTaggingHistos() ;
613
614 //Remove old events
615 fCurrentMixedList->AddFirst(fPHOSEvent);
616 fPHOSEvent=0x0 ;
617 if(fCurrentMixedList->GetSize() > 10){
618 TClonesArray *tmp = static_cast <TClonesArray*> (fCurrentMixedList->Last());
619 fCurrentMixedList->Remove(tmp);
620 delete tmp;
621 }
622
623 PostData(1, fOutputContainer);
624
625}
626//________________________________________________
627void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
628
629 //MC info about this particle
630 if(!fStack)
631 return ;
632
633 const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
634 AliAODEvent* event = (AliAODEvent*)InputEvent();
635
636 //Fill spectra of primary particles
637 char partName[10] ;
638 for(Int_t i=0;i<fStack->GetNtrack();i++){
639 TParticle* particle = fStack->Particle(i);
640 //Primary particle
641 Double_t r=particle->R() ;
642
643 switch(particle->GetPdgCode()){
644 case 111: snprintf(partName,10,"pi0") ;
645 break ;
646 case 221: snprintf(partName,10,"eta") ;
647 break ;
648 case 223: snprintf(partName,10,"omega") ;
649 break ;
650 case 22: snprintf(partName,10,"gamma") ;
651 break ;
652 default: continue ;
653 }
654
655 Double_t pt = particle->Pt() ;
656 //Distribution over vertex
657 FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
658
659 if(r >rcut)
660 continue ;
661
662 //Total number of pi0 with creation radius <1 cm
663 Double_t w = PrimaryParticleWeight(particle) ;
664 FillHistogram(Form("hMC_all_%s",partName),pt,w) ;
665 if(TMath::Abs(particle->Y())<0.15){
666 FillHistogram(Form("hMC_unitEta_%s",partName),pt,w) ;
667 }
668
669 FillHistogram(Form("hMC_rap_%s",partName),particle->Y(),w) ;
670
671 Double_t phi=particle->Phi() ;
672 while(phi<0.)phi+=TMath::TwoPi() ;
673 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
674 FillHistogram(Form("hMC_phi_%s",partName),phi,w) ;
675 }
676
677
678 //Clussify reconstructed clusters
679 //First - photons (from vertex) and contaminations
680 //Second - photons from different sources
681 //Third - photons from pi0s - missed for different reasons
682
683 const Int_t n=fPHOSEvent->GetEntriesFast() ;
684 for(Int_t i=0;i<n;i++){
685 AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
686 Int_t iprim=p->GetCaloLabel(0) ; //Particle entered PHOS
687 if(iprim<0){ //No label!
688 FillHistogram("hMCRecNoLabel",p->Pt());
689 continue ;
690 }
691
692 TParticle * prim = fStack->Particle(iprim) ;
693
694 //Look what particle left vertex
695 Int_t iparent=p->GetCaloLabel(1); //Particle left vertex
696 TParticle * parent =fStack->Particle(iparent);
697
698 Int_t grandpaPDG=-1 ;
699 if(parent){
700 grandpaPDG=parent->GetPdgCode() ;
701 }
702 switch(grandpaPDG){
703 case -1: //no primary
704 FillPIDHistograms("hMCRecNoPrim",p);
705 break ;
706 case 22:
707 FillPIDHistograms("hMCRecGamma",p);
708 break ;
709 case 11:
710 case -11: //electron/positron conversion
711 FillPIDHistograms("hMCRecPhotConv",p); //Reconstructed with photon from conversion primary
712 FillHistogram("hMCConversionRadius",prim->R());
713 break ;
714 case 111: //Pi0 decay //Primary decay photon (as in MC)
715 FillPIDHistograms("hMCRecPi0",p);
716 //Strange, vertex of pi0?
717 FillHistogram("hMCRecPi0Vtx",p->Pt(),parent->R());
718 break ;
719 case 221: //eta decay
720 FillPIDHistograms("hMCRecEta",p);
721 //Strange, vertex of eta?
722 FillHistogram("hMCRecEtaVtx",p->Pt(),parent->R());
723 break ;
724 case 223: //omega meson decay
725 FillPIDHistograms("hMCRecOmega",p);
726 //Strange, vertex of omega?
727 FillHistogram("hMCRecOmegaVtx",p->Pt(),parent->R());
728 break ;
729 case 331: //eta' decay
730 FillPIDHistograms("hMCRecEtapr",p);
731 //Strange, vertex of eta'?
732 FillHistogram("hMCRecEtaprVtx",p->Pt(),parent->R());
733 break ;
734 case -2212:
735 FillPIDHistograms("hMCRecPbar",p); //Reconstructed with photon from antibaryon annihilation
736 break ;
737 case -2112: //antineutron & antiproton
738 FillPIDHistograms("hMCRecNbar",p); //Reconstructed with photon from antibaryon annihilation
739 break ;
740 case 211:
741 case -211:
742 FillPIDHistograms("hMCRecPipm",p); //Reconstructed with photon from antibaryon annihilation
743 break ;
744 case 2112:
745 FillPIDHistograms("hMCRecN",p); //Reconstructed with photon from antibaryon annihilation
746 break ;
747 case 2212:
748 FillPIDHistograms("hMCRecP",p); //Reconstructed with photon from antibaryon annihilation
749 break ;
750 case 321:
751 case -321:
752 FillPIDHistograms("hMCRecKpm",p); //Reconstructed with photon from conversion primary
753 break ;
754 case 130:
755 FillHistogram("hMCRecK0sVtx",p->Pt(),prim->R());
756 FillPIDHistograms("hMCRecK0s",p); //Reconstructed with photon from conversion primary
757 break ;
758 case 310:
759 FillHistogram("hMCRecK0lVtx",p->Pt(),prim->R());
760 FillPIDHistograms("hMCRecK0l",p); //Reconstructed with photon from conversion primary
761 break ;
762 default:
763 if(parent->GetPDG()->Charge()!=0)
764 FillPIDHistograms("hMCRecUnknownCh",p);
765 else
766 FillPIDHistograms("hMCRecUnknownNeu",p);
767 printf("Unknown PDG code: %d \n",grandpaPDG) ;
768 }
769//Put here filling histograms vs PID...
770
771
772 //Now classify decay photons
773 if(grandpaPDG==22){
774 //Separate gammas from pi0s, eta, omega, eta', e+-, other, no primary
775 Int_t igrandpa = -1 ;
776 if(parent)
777 igrandpa=parent->GetFirstMother();
778 Int_t decayPDG=-1 ;
779 TParticle * decay=0x0 ;
780 if(igrandpa>-1){
781 decay=fStack->Particle(igrandpa) ;
782 decayPDG=decay->GetPdgCode() ;
783 }
784 switch(decayPDG){
785 case -1 :// Direct gamma?
786 FillPIDHistograms("hMCRecGammaDir",p);
787 break ;
788 case 111 :// pi0 decays
789 FillPIDHistograms("hMCRecGammaPi0",p);
790 break ;
791 case 221 :// eta decays
792 FillPIDHistograms("hMCRecGammaEta",p);
793 break ;
794 case 223 :// omega decays
795 FillPIDHistograms("hMCRecGammaOmega",p);
796 break ;
797 default:
798 FillPIDHistograms("hMCRecGammaOther",p);
799 }
800
801
802 //Classify photons from pi0 decay
803 if(decayPDG==111){
804 //Dalitz decay
805 if(decay->GetNDaughters()>2){
806 FillPIDHistograms("hMCRecGammaPi0Dalitz",p);
807 continue ;
808 }
809 Int_t ipartner = decay->GetFirstDaughter();
810 if(ipartner==iparent){ //look other
811 ipartner = decay->GetLastDaughter();
812 }
813 //No partner in Stack
814 if(ipartner==-1){
815 FillPIDHistograms("hMCRecGammaPi0NoStack",p) ;
816 }
817 else{
818 //
819 FillPIDHistograms("hMCGammaPi02Gamma",p) ;
820 TParticle * partner = fStack->Particle(ipartner);
821 //Check if partner is registered and made correct mass
822 Double_t mPrim=(parent->Energy()+partner->Energy())*(parent->Energy()+partner->Energy())-
823 (parent->Px()+partner->Px())*(parent->Px()+partner->Px()) -
824 (parent->Py()+partner->Py())*(parent->Py()+partner->Py()) -
825 (parent->Pz()+partner->Pz())*(parent->Pz()+partner->Pz()) ;
826 if(mPrim>0.)mPrim=TMath::Sqrt(mPrim) ;
b5f58f9c 827 FillHistogram(Form("hMCGammaPi0PrimMgg_cent%d",fCentBin),mPrim,p->Pt()) ;
213dd723 828
829 //find corresponding PWG4particle
830 Bool_t isPartnerRec=kFALSE ;
831 Bool_t isTagged=kFALSE ;
832 for(Int_t j=0;j<n;j++){
833 if(j==i)
834 continue ;
835 AliAODPWG4Particle *p2 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(j));
836 Double_t invMass = p->GetPairMass(p2);
837 if(p2->GetCaloLabel(1)==ipartner){ //partner
838 isPartnerRec=kTRUE ;
839 FillHistogram("hMCGammaPi0RecMgg",invMass,p->Pt()) ;
840 FillPIDHistograms("hMCGammaPi0Rec",p) ;
841
842 //estimate proportion passed PI0 cut
843 if(IsInPi0Band(invMass,p->Pt(),2)){ //Check type!!!!!!!!!!!!!!
844 FillPIDHistograms("hMCGammaPi0TrueTagged",p) ;
845 if(isTagged)
846 FillPIDHistograms("hMCGammaPi0MultyTagged",p) ;
847 isTagged=kTRUE ;
848 }
849 }
850 else{
851 if(IsInPi0Band(invMass,p->Pt(),2)){//Check type!!!!!!!!!!!!!!
852 FillPIDHistograms("hMCGammaPi0FakeTagged",p) ;
853 if(isTagged)
854 FillPIDHistograms("hMCGammaPi0MultyTagged",p) ;
855 isTagged=kTRUE ;
856 }
857 }
858 }
859 if(isTagged)
860 FillPIDHistograms("hMCGammaPi0Tagged",p) ;
861
862 //If no PWG4particle trace the reason
863 if(!isPartnerRec){
864 //conversion
865 if(partner->GetPdgCode()==22){
866 Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
867
868 //did not hit PHOS
869 Int_t modulenum ;
870 Double_t ztmp=0.,xtmp=0. ;
871 Bool_t impact=fPHOSgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp) ;
872 //Check Bad Map
873 Int_t fidArea=0 ;
874 TVector3 partnerGlobaPos ;
875 if(impact){
876 fPHOSgeom->Local2Global(modulenum,xtmp, ztmp,partnerGlobaPos) ;
877 Float_t pos[3] ;
878 partnerGlobaPos.GetXYZ(pos);
879 fidArea=GetFiducialArea(pos) ;
880 }
881
882 if(!impact || (fidArea==0) ){ //this photon cannot hit PHOS
883 FillPIDHistograms("hMCGammaPi0MisGeo",p);
884 Int_t iFidArea = p->GetFiducialArea();
885 if(iFidArea>0){
886 FillPIDHistograms("hMCGammaPi0MisGeoFA1",p) ; //Spectrum of tagged with missed partner
887 if(iFidArea>1){
888 FillPIDHistograms("hMCGammaPi0MisGeoFA2",p) ; //Spectrum of tagged with missed partner
889 if(iFidArea>2){
890 FillPIDHistograms("hMCGammaPi0MisGeoFA3",p) ; //Spectrum of tagged with missed partner
891 }
892 }
893 }
894 isPartnerLost=kTRUE;
895 }
896
897 //this photon is converted before it is registered
898 if(!isPartnerLost){
899 if(partner->GetNDaughters()>0 &&
900 fStack->Particle(partner->GetFirstDaughter())->R()<450.){
901 FillPIDHistograms("hMCGammaPi0MisConv",p);
902 FillHistogram("hMCGammaPi0MisConvR",p->Pt(),fStack->Particle(partner->GetFirstDaughter())->R()) ; //Spectrum of tagged with missed partner
903 isPartnerLost=kTRUE;
904 }
905 }
906
907 // too soft
908 if(!isPartnerLost &&
909 partner->Energy()<0.3){ //energy is not enough to be registered by PHOS
910 FillPIDHistograms("hMCGammaPi0MisEmin",p);
911 isPartnerLost=kTRUE;
912 }
913 //Other reason
914 if(!isPartnerLost){
915 FillHistogram("hMCGammaPi0MisPartner",partner->Pt());
916 FillHistogram("hMCGammaPi0MisPartnerEtaPhi",partner->Eta(),partner->Phi());
917 //May be overlap?
918 Bool_t fakePrimary=kFALSE ;
919 Int_t multClust = event->GetNumberOfCaloClusters();
920 for (Int_t iclu=0; iclu<multClust; iclu++) {
921 AliAODCaloCluster * clu = event->GetCaloCluster(iclu);
922 Float_t pos[3] ;
923 clu->GetPosition(pos) ;
924 TVector3 global1(pos) ;
925 Double_t d=(partnerGlobaPos-global1).Mag();
926 if( d<4 ){//same cluster
927 //check is partner is in list?
928 Int_t nlab=clu->GetNLabels() ;
929 for(Int_t ilab=0; ilab<nlab; ilab++){
930 Int_t labelA = clu->GetLabelAt(ilab) ;
931 while(labelA>-1){
932 if(labelA==ipartner){
933 if(clu->E()>0.3)
934 FillPIDHistograms("hMCGammaPi0RecSoft",p);
935 else
936 if(clu->GetNCells()<3)
937 FillPIDHistograms("hMCGammaPi0RecCells",p);
938 else
939 FillPIDHistograms("hMCGammaPi0MisFoundPrim",p);
940 break ;
941 }
942 labelA=fStack->Particle(labelA)->GetFirstMother() ;
943 }
944 }
945 fakePrimary=kTRUE ;
946 break ;
947 }
948 }//end of loop over clusters
949 //No cluster in this region
950 if(fakePrimary)
951 FillPIDHistograms("hMCGammaPi0MisFakePrim",p);
952 else
953 FillPIDHistograms("hMCGammaPi0MisOther",p);
954
955 }
956 }//Partner - photon
957 else{//partner not photon
958 FillPIDHistograms("hMCGammaPi0MisNPhot",p);
959 }
960 }
961 }
962 }
963 }
964 } //PHOS clusters
965
966
967 //Fill histograms with all clusters fake tagging
968 for(Int_t i=0;i<n-1;i++){
969 AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
970
971 Int_t ipartner=-999 ; // this will be partner if this cluster from pi0 decay
972 //Look what particle left vertex
973 Int_t iparent=p->GetCaloLabel(1); //Particle left vertex
974 TParticle * parent =fStack->Particle(iparent);
975 if(parent){
976 Int_t grandpaPDG=parent->GetPdgCode() ;
977 if(grandpaPDG==22){
978 Int_t igrandpa =parent->GetFirstMother();
979 if(igrandpa>-1){
980 TParticle * decay=fStack->Particle(igrandpa) ;
981 Int_t decayPDG=decay->GetPdgCode() ;
982 if(decayPDG==111){
983 //Dalitz decay
984 if(decay->GetNDaughters()==2){
985 ipartner = decay->GetFirstDaughter();
986 if(ipartner==iparent){ //look other
987 ipartner = decay->GetLastDaughter();
988 }
989 if(ipartner==-1) //no partner
990 ipartner=-999 ;
991 }
992 }
993 }
994 }
995 }
996 for(Int_t j=i+1;j<n;j++){
997 AliAODPWG4Particle *p2 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(j));
998 if(p2->GetCaloLabel(1)!=ipartner){ //partner
999 Double_t invMass = p->GetPairMass(p2);
1000 if(IsInPi0Band(invMass,p->Pt(),2)){
1001 FillPIDHistograms("hMCAllFakeTagged",p) ;
1002 }
1003 if(IsInPi0Band(invMass,p2->Pt(),2)){
1004 FillPIDHistograms("hMCAllFakeTagged",p2) ;
1005 }
1006 }
1007 }
1008 }
1009
1010}
1011//________________________________________________
1012void AliAnalysisTaskTaggedPhotons::FillTaggingHistos(){
1013 //Fill all necessary histograms
1014
1015
1016 //First fill invariant mass distrubutions and mark tagged photons
1017 //Invariant Mass analysis
1018 const Int_t n=fPHOSEvent->GetEntriesFast() ;
1019 for(Int_t i=0;i<n-1;i++){
1020 AliAODPWG4Particle *p1 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
1021 for(Int_t j = i+1 ; j < n ; j++) {
1022 AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(j));
1023
1024 Double_t invMass = p1->GetPairMass(p2);
1025
1026 if(p2->E()>0.1){
1027 FillPIDHistograms("hInvM_Re_Emin1",p1,invMass) ;
1028 if(p2->E()>0.2){
1029 FillPIDHistograms("hInvM_Re_Emin2",p1,invMass) ;
1030 if(p2->E()>0.3){
1031 FillPIDHistograms("hInvM_Re_Emin3",p1,invMass) ;
1032 }
1033 }
1034 }
1035
1036 if(p1->E()>0.1){
1037 FillPIDHistograms("hInvM_Re_Emin1",p2,invMass) ;
1038 if(p1->E()>0.2){
1039 FillPIDHistograms("hInvM_Re_Emin2",p2,invMass) ;
1040 if(p1->E()>0.3){
1041 FillPIDHistograms("hInvM_Re_Emin3",p2,invMass) ;
1042 }
1043 }
1044 }
1045 if(TestPID(3, p2)){
1046 if(p2->E()>0.1){
1047 FillPIDHistograms("hInvM_Re_Emin1_Strict",p1,invMass) ;
1048 if(p2->E()>0.2){
1049 FillPIDHistograms("hInvM_Re_Emin2_Strict",p1,invMass) ;
1050 if(p2->E()>0.3){
1051 FillPIDHistograms("hInvM_Re_Emin3_Strict",p1,invMass) ;
1052 }
1053 }
1054 }
1055 }
1056 if(TestPID(3, p1)){
1057 if(p1->E()>0.1){
1058 FillPIDHistograms("hInvM_Re_Emin1_Strict",p2,invMass) ;
1059 if(p1->E()>0.1){
1060 FillPIDHistograms("hInvM_Re_Emin2_Strict",p2,invMass) ;
1061 if(p1->E()>0.3){
1062 FillPIDHistograms("hInvM_Re_Emin3_Strict",p2,invMass) ;
1063 }
1064 }
1065 }
1066 }
1067 if(IsSamePi0(p1,p2)){
1068 FillPIDHistograms("hMC_InvM_Re",p1,invMass) ;
1069 FillPIDHistograms("hMC_InvM_Re",p2,invMass) ;
1070 if(TestPID(3, p2)){
1071 FillPIDHistograms("hMC_InvM_Re_Strict",p1,invMass) ;
1072 }
1073 if(TestPID(3, p1)){
1074 FillPIDHistograms("hMC_InvM_Re_Strict",p2,invMass) ;
1075 }
1076 }
1077
1078 //Tagging: 1,2,3 sigma
1079 //Emin=100,200,300 Mev
1080 //IsInPi0Band(mass, Ptphoton, type Emin cut
1081 Int_t tag1=0 ;
1082 for(Int_t eminType=0; eminType<3; eminType++){
1083 if(p2->E()>0.1*(eminType+1)){
1084 //Set corresponding bit
1085 Double_t nsigma = IsInPi0Band(invMass,p1->Pt(),eminType) ; //in band with n sigmas
1086 for(Int_t isigma=0; isigma<3; isigma++){
1087 if(nsigma<1+isigma){
1088 tag1|= (1 << (3*eminType+isigma)) ;
1089 if(TestPID(3, p2))
1090 tag1|= (1 << (3*eminType+isigma+9)) ;
1091 }
1092 }
1093 }
1094 }
1095 p1->SetTag(tag1) ;
1096 Int_t tag2=0 ;
1097 for(Int_t eminType=0; eminType<3; eminType++){
1098 if(p1->E()>0.1*(eminType+1)){
1099 //Set corresponding bit
1100 Double_t nsigma = IsInPi0Band(invMass,p2->Pt(),eminType) ; //in band with n sigmas
1101 for(Int_t isigma=0; isigma<3; isigma++){
1102 if(nsigma<1+isigma){
1103 tag2|= (1 << (3*eminType+isigma)) ;
1104 if(TestPID(3, p2))
1105 tag2|= (1 << (3*eminType+isigma+9)) ;
1106 }
1107 }
1108 }
1109 }
1110 p2->SetTag(tag2) ;
1111
1112 if(tag1 & (1<<7)){ //2 sigma, Emin=0.3: default tagging
1113 if(p1->IsTagged()){//Multiple tagging
1114 FillHistogram(Form("hTaggedMult_cent%d",fCentBin),p1->Pt());
1115 }
1116 p1->SetTagged(kTRUE) ;
1117 }
1118 if(tag2 & (1<<7)){ //2 sigma, Emin=0.3: default tagging
1119 if(p2->IsTagged()){//Multiple tagging
1120 FillHistogram(Form("hTaggedMult_cent%d",fCentBin),p2->Pt());
1121 }
1122 p2->SetTagged(kTRUE) ;
1123 }
1124 }
1125 }
1126
1127
1128
1129 //Single particle histgams
1130 for(Int_t i=0;i<n;i++){
1131 AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
1132
1133 Int_t isolation = p->GetBtag();
1134
1135 //Inclusive spectra
1136 FillPIDHistograms("hPhot",p) ;
1137
1138 if(p->DistToBad()>1){
1139 FillPIDHistograms("hPhot_Dist2",p) ;
1140 if(p->DistToBad()>2){
1141 FillPIDHistograms("hPhot_Dist3",p) ;
1142 }
1143 }
1144
1145
1146 for(Int_t kind=1; kind<33; kind*=2){
1147 if((isolation&kind)){
1148 FillPIDHistograms(Form("hPhot_Isolation%d",kind),p) ;
1149 }
1150 }
1151
1152 Int_t iFidArea = p->GetFiducialArea();
1153 if(iFidArea>0){
1154 FillPIDHistograms("hPhot_Area1",p) ;
1155 for(Int_t kind=1; kind<33; kind*=2){
1156 if((isolation&kind)){
1157 FillPIDHistograms(Form("hPhot_Isolation%d_Area1",kind),p) ;
1158 }
1159 }
1160
1161 //Fill taggings with
1162 //3 Emin cuts
1163 //Default Emin, 1,2,3 sigmas
1164 //strict and loose PID cut on partner
1165 Int_t tag=p->GetTag() ;
1166 for(Int_t ibit=0; ibit<18; ibit++){
1167 if((tag & (1<<ibit))==0){
1168 FillPIDHistograms(Form("hPhot_nTagged%d_Area1",ibit),p) ;
1169// for(Int_t kind=1; kind<33; kind*=2){
1170// if((isolation&kind)){
1171// FillPIDHistograms(Form("hPhot_nStrictTagged_Isolation%d_Area1",kind),p) ;
1172// }
1173// }
1174 }
1175 }
1176
1177 if(iFidArea>1){
1178 FillPIDHistograms("hPhot_Area2",p) ;
1179 for(Int_t ibit=0; ibit<18; ibit++){
1180 if((tag & (1<<ibit))==0){
1181 FillPIDHistograms(Form("hPhot_nTagged%d_Area2",ibit),p) ;
1182 }
1183 }
1184 if(iFidArea>2){
1185 FillPIDHistograms("hPhot_Area3",p) ;
1186 for(Int_t ibit=0; ibit<18; ibit++){
1187 if((tag & (1<<ibit))==0){
1188 FillPIDHistograms(Form("hPhot_nTagged%d_Area3",ibit),p) ;
1189 }
1190 }
1191 }
1192 }
1193 }
1194 }
1195
1196 //Fill Mixed InvMass distributions:
1197 TIter nextEv(fCurrentMixedList) ;
1198 for(Int_t i=0;i<n;i++){
1199 AliAODPWG4Particle *p1 = static_cast<AliAODPWG4Particle*>(fPHOSEvent->At(i));
1200 while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
1201 Int_t nPhotons2 = event2->GetEntriesFast() ;
1202 for(Int_t j=0; j < nPhotons2 ; j++){
1203 AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;
1204 Double_t invMass = p1->GetPairMass(p2);
1205
1206 if(p2->E()>0.1){
1207 FillPIDHistograms("hInvM_Mi_Emin1",p1,invMass) ;
1208 if(p2->E()>0.2){
1209 FillPIDHistograms("hInvM_Mi_Emin2",p1,invMass) ;
1210 if(p2->E()>0.3){
1211 FillPIDHistograms("hInvM_Mi_Emin3",p1,invMass) ;
1212 }
1213 }
1214 }
1215 if(TestPID(3, p2)){
1216 if(p2->E()>0.1){
1217 FillPIDHistograms("hInvM_Mi_Emin1_Strict",p1,invMass) ;
1218 if(p2->E()>0.2){
1219 FillPIDHistograms("hInvM_Mi_Emin2_Strict",p1,invMass) ;
1220 if(p2->E()>0.3){
1221 FillPIDHistograms("hInvM_Mi_Emin3_Strict",p1,invMass) ;
1222 }
1223 }
1224 }
1225 }
1226
1227 if(p1->E()>0.1){
1228 FillPIDHistograms("hInvM_Mi_Emin1",p2,invMass) ;
1229 if(p1->E()>0.2){
1230 FillPIDHistograms("hInvM_Mi_Emin2",p2,invMass) ;
1231 if(p1->E()>0.3){
1232 FillPIDHistograms("hInvM_Mi_Emin3",p2,invMass) ;
1233 }
1234 }
1235 }
1236 if(TestPID(3, p1)){
1237 if(p1->E()>0.1){
1238 FillPIDHistograms("hInvM_Mi_Emin1_Strict",p2,invMass) ;
1239 if(p1->E()>0.2){
1240 FillPIDHistograms("hInvM_Mi_Emin2_Strict",p2,invMass) ;
1241 if(p1->E()>0.3){
1242 FillPIDHistograms("hInvM_Mi_Emin3_Strict",p2,invMass) ;
1243 }
1244 }
1245 }
1246 }
1247 }
1248 }
1249 }
1250
1251
1252}
1253
1254//______________________________________________________________________________
1255void AliAnalysisTaskTaggedPhotons::Init()
1256{
1257 // Intialisation of parameters
1258}
1259
1260//______________________________________________________________________________
1261void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
1262{
1263 // Processing when the event loop is ended
1264 if (fDebug > 1) Printf("Terminate()");
1265}
1266//______________________________________________________________________________
1267Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt, Int_t /*type*/)const
1268{
1269 //Parameterization of the pi0 peak region
1270 Double_t mpi0mean = 1.33259e-01 - 2.910e-03 * TMath::Exp(-pt/3.616) ;
1271// Double_t mpi0mean = 0.136 ; //fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;
1272
1273 Double_t mpi0sigma=TMath::Sqrt(4.17927e-03*4.17927e-03+2.81581e-03*2.81581e-03/pt+3.59218e-04*3.59218e-04*pt*pt) ;
1274// Double_t mpi0sigma= 0.005; //TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);
1275
1276 return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
1277}
1278//______________________________________________________________________________
1279Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{
1280 //Looks through parents and finds if there was commont pi0 among ancestors
1281
1282 if(!fStack)
1283 return kFALSE ; //can not say anything
1284
1285 Int_t prim1 = p1->GetCaloLabel(0);
1286 while(prim1!=-1){
1287 Int_t prim2 = p2->GetCaloLabel(0);
1288 while(prim2!=-1){
1289 if(prim1==prim2){
1290 if(fStack->Particle(prim1)->GetPdgCode()==111)
1291 return kTRUE ;
1292 else
1293 return kFALSE ;
1294 }
1295 prim2=fStack->Particle(prim2)->GetFirstMother() ;
1296 }
1297 prim1=fStack->Particle(prim1)->GetFirstMother() ;
1298 }
1299 return kFALSE ;
1300}
1301//______________________________________________________________________________
1302Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(const Float_t * position)const{
1303 //calculates in which kind of fiducial area photon hit
1304
1305 TVector3 global1(position) ;
1306 Int_t relId[4] ;
1307 fPHOSgeom->GlobalPos2RelId(global1,relId) ;
1308 Int_t mod = relId[0] ;
1309 Int_t cellX = relId[2];
1310 Int_t cellZ = relId[3] ;
1311
1312 if(fPHOSBadMap[mod] && fPHOSBadMap[mod]->GetBinContent(cellX,cellZ)>0)
1313 return 0 ;
1314
1315 //New we are in good channel,
1316 //calculate distance to the closest group of bad channels
1317 const Int_t cut1=10;
1318 const Int_t cut2=15;
1319//For 3-mod configuration
1320// if((mod==3 && cellX<=cut1) || (mod==1 && cellX>=65-cut1) || cellZ<=cut1 || cellZ>=57-cut1)
1321//For 1+3 configuraion
1322 if( cellX<=cut1 || cellX>=65-cut1 || cellZ<=cut1 || cellZ>=57-cut1)
1323 return 1;
1324// //and from large dead area
1325// if(fPHOSBadMap[mod]->Integral(cellX-cut1,cellX+cut1,cellZ-cut1,cellZ+cut1)>0.2*cut1*cut1)
1326// return 1 ;
1327//Full configuration
1328// if((mod==3 && cellX<=cut2) || (mod==1 && cellX>=65-cut2) || cellZ<=cut2 || cellZ>=57-cut2)
1329//1+3 configuration
1330 if( cellX<=cut2 || cellX>=65-cut2 || cellZ<=cut2 || cellZ>=57-cut2)
1331 return 2;
1332// if(fPHOSBadMap[mod]->Integral(cellX-cut2,cellX+cut2,cellZ-cut2,cellZ+cut2)>0.8*cut2*cut2)
1333// return 2 ;
1334 //Very good channel
1335 return 3 ;
1336
1337}
1338//______________________________________________________________________________^M
1339void AliAnalysisTaskTaggedPhotons::InitGeometry(){
1340 //Rotation matrixes are set with Tender
1341
1342 if(fPHOSgeom) return ;
1343 else
1344 fPHOSgeom = AliPHOSGeometry::GetInstance() ;
1345
1346}
1347//_____________________________________________________________________________
1348void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x)const{
1349 //FillHistogram
1350 TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
1351 if(tmpI){
1352 tmpI->Fill(x) ;
1353 return ;
1354 }
1355 TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
1356 if(tmpF){
1357 tmpF->Fill(x) ;
1358 return ;
1359 }
1360 TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
1361 if(tmpD){
1362 tmpD->Fill(x) ;
1363 return ;
1364 }
1365 AliInfo(Form("can not find histogram <%s> ",key)) ;
1366}
1367//_____________________________________________________________________________
1368void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x,Double_t y)const{
1369 //FillHistogram
1370 TObject * tmp = fOutputContainer->FindObject(key) ;
1371 if(!tmp){
1372 AliInfo(Form("can not find histogram <%s> ",key)) ;
1373 return ;
1374 }
1375 if(tmp->IsA() == TClass::GetClass("TH1F")){
1376 ((TH1F*)tmp)->Fill(x,y) ;
1377 return ;
1378 }
1379 if(tmp->IsA() == TClass::GetClass("TH2F")){
1380 ((TH2F*)tmp)->Fill(x,y) ;
1381 return ;
1382 }
1383 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1384}
1385
1386//_____________________________________________________________________________
1387void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1388 //Fills 1D histograms with key
1389 TObject * tmp = fOutputContainer->FindObject(key) ;
1390 if(!tmp){
1391 AliInfo(Form("can not find histogram <%s> ",key)) ;
1392 return ;
1393 }
1394 if(tmp->IsA() == TClass::GetClass("TH2F")){
1395 ((TH2F*)tmp)->Fill(x,y,z) ;
1396 return ;
1397 }
1398 if(tmp->IsA() == TClass::GetClass("TH3F")){
1399 ((TH3F*)tmp)->Fill(x,y,z) ;
1400 return ;
1401 }
1402}
1403//_____________________________________________________________________________
1404void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliAODPWG4Particle * p) const{
1405
1406 FillHistogram(Form("%s_All_cent%d",name,fCentBin),p->Pt()) ;
1407 if(p->GetDispBit())
1408 FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),p->Pt()) ;
1409 if(p->GetChargedBit())
1410 FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),p->Pt()) ;
1411 if(p->GetDispBit() && p->GetChargedBit())
1412 FillHistogram(Form("%s_Both_cent%d",name,fCentBin),p->Pt()) ;
1413
1414}
1415//_____________________________________________________________________________
1416void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliAODPWG4Particle * p,Double_t x) const{
1417
1418 FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,p->Pt()) ;
1419 if(p->GetDispBit())
1420 FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,p->Pt()) ;
1421 if(p->GetChargedBit())
1422 FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,p->Pt()) ;
1423 if(p->GetDispBit() && p->GetChargedBit())
1424 FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,p->Pt()) ;
1425
1426}
1427//_____________________________________________________________________________
1428Bool_t AliAnalysisTaskTaggedPhotons::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1429
1430 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1431 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1432 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1433 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1434 Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1435 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1436 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1437 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1438 return (R2<2.5*2.5) ;
1439
1440}
1441//_________________________________________________________________________________
1442Int_t AliAnalysisTaskTaggedPhotons::EvalIsolation(TLorentzVector * ph){
1443
1444 // Check if this particle is isolated.
1445 //We use several cone radii and epsilons.
1446 //As well we look at charged particles and EMCAL ones
1447
1448 const Double_t coneR1=0.3 ;
1449 const Double_t coneR2=0.4 ;
1450 const Double_t coneR3=0.5 ;
1451
1452 const Double_t epsilon1=0.1 ;
1453 const Double_t epsilon2=0.05 ;
1454
1455 if(!ph) return 0 ;
1456
1457 //Sum of energies in cones, tracks and clusters in EMCAL
1458 Double_t eCone1 = 0;
1459 Double_t eCone2 = 0;
1460 Double_t eCone3 = 0;
1461
1462 Double_t phiTrig = ph->Phi();
1463 Double_t etaTrig = ph->Eta();
1464
1465 Int_t n=fTrackEvent->GetEntriesFast() ;
1466 for(Int_t itr=0; itr<n; itr++){
1467 AliAODPWG4Particle * track = (AliAODPWG4Particle*)fTrackEvent->At(itr) ;
1468
1469 Double_t deleta = etaTrig - track->Eta() ;
1470 Double_t delphi = phiTrig - track->Phi() ;
1471 while(delphi<-TMath::Pi()) delphi+=TMath::TwoPi() ;
1472 while(delphi>TMath::Pi()) delphi-=TMath::TwoPi() ;
1473
1474 Double_t dr = TMath::Sqrt(deleta * deleta + delphi * delphi);
1475
1476 if(dr<coneR3){
1477 eCone3+=track->Pt() ;
1478 if(dr<coneR2){
1479 eCone2+=track->Pt() ;
1480 if(dr<coneR1){
1481 eCone1+=track->Pt() ;
1482 }
1483 }
1484 }
1485 }
1486
1487 //Fill QA histgams
1488 Double_t ptTrig=ph->Pt() ;
1489 FillHistogram(Form("QA_Cone1_Tracks_cent%d",fCentBin),ptTrig,eCone1) ;
1490 FillHistogram(Form("QA_Cone2_Tracks_cent%d",fCentBin),ptTrig,eCone2) ;
1491 FillHistogram(Form("QA_Cone3_Tracks_cent%d",fCentBin),ptTrig,eCone3) ;
1492
1493
1494 //Fill Bits
1495 Int_t iCone1E1 = (epsilon1*ptTrig > eCone1) ;
1496 Int_t iCone2E1 = (epsilon1*ptTrig > eCone2) ;
1497 Int_t iCone3E1 = (epsilon1*ptTrig > eCone3) ;
1498
1499 Int_t iCone1E2 = (epsilon2*ptTrig > eCone1) ;
1500 Int_t iCone2E2 = (epsilon2*ptTrig > eCone2) ;
1501 Int_t iCone3E2 = (epsilon2*ptTrig > eCone3) ;
1502
1503
1504 Int_t isolation= iCone1E1+ 2*iCone2E1 +4*iCone3E1+
1505 8*iCone1E2 +16*iCone2E2 +32*iCone3E2 ;
1506 return isolation ;
1507}
1508//_________________________________________________________________________________
1509Bool_t AliAnalysisTaskTaggedPhotons::TestPID(Int_t iPID, AliAODPWG4Particle* part){
1510// //Checks PID of particle
1511
1512 if(!part) return kFALSE ;
1513 if(iPID==0) return kTRUE ;
1514 if(iPID==1) return part->GetDispBit() ;
1515 if(iPID==2) return part->GetChargedBit() ;
1516 if(iPID==3) return part->GetDispBit() && part->GetChargedBit() ;
1517 return kFALSE ;
1518
1519}
1520//_________________________________________________________________________________
1521Double_t AliAnalysisTaskTaggedPhotons::PrimaryParticleWeight(TParticle * /*particle*/){
1522 return 1;
1523
1524}
1525//___________________________________________________________________________
b5f58f9c 1526Int_t AliAnalysisTaskTaggedPhotons::FindPrimary(AliVCluster*clu, Bool_t&sure){
213dd723 1527 //Finds primary and estimates if it unique one?
1528 //First check can it be photon/electron
1529 const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
1530 Int_t n=clu->GetNLabels() ;
b5f58f9c 1531 FillHistogram(Form("LabelsNPrim_cent%d",fCentBin),clu->E(),n) ;
213dd723 1532 if(n==1){
1533 sure=kTRUE;
1534 return clu->GetLabelAt(0) ;
1535 }
1536 //Number of clusters with one or more photons
1537 Bool_t hasGamma=kFALSE ;
1538 Double_t eMax=0. ;
1539 for(Int_t i=0; i<n; i++){
1540 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
1541 Int_t pdg = p->GetPdgCode() ;
1542 if(pdg==22){
1543 hasGamma=kTRUE ;
1544 if(p->Energy()>eMax){
1545 eMax=p->Energy();
1546 }
1547 }
1548 }
1549 if(hasGamma){
b5f58f9c 1550 FillHistogram(Form("LabelsGamma_cent%d",fCentBin),clu->E()) ;
1551 FillHistogram(Form("LabelsGammaE_cent%d",fCentBin),clu->E(),eMax/clu->E()) ;
213dd723 1552 }
1553
1554 for(Int_t i=0; i<n; i++){
1555 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
1556 Int_t pdg = p->GetPdgCode() ;
1557 if(pdg==22 || pdg==11 || pdg == -11){
1558 if(p->Energy()>emFraction*clu->E()){
1559 sure=kTRUE ;
1560 return clu->GetLabelAt(i);
1561 }
1562 }
1563 }
1564
1565 Double_t* Ekin= new Double_t[n] ;
1566 for(Int_t i=0; i<n; i++){
1567 TParticle* p= fStack->Particle(clu->GetLabelAt(i)) ;
1568 Ekin[i]=p->P() ; // estimate of kinetic energy
1569 if(p->GetPdgCode()==-2212 || p->GetPdgCode()==-2112){
1570 Ekin[i]+=1.8 ; //due to annihilation
1571 }
1572 }
1573 Int_t iMax=0;
1574 Double_t eSubMax=0. ;
1575 eMax=0.;
1576 for(Int_t i=0; i<n; i++){
1577 if(Ekin[i]>eMax){
1578 eSubMax=eMax;
1579 eMax=Ekin[i];
1580 iMax=i;
1581 }
1582 }
1583 if(eSubMax>0.8*eMax)//not obvious primary
1584 sure=kFALSE;
1585 else
1586 sure=kTRUE;
1587 delete[] Ekin;
1588 return clu->GetLabelAt(iMax) ;
1589}
1590
1591