]> 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
Weights for MC simulations added
[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"
0b9a6a23 44#include "AliCaloPhoton.h"
92071be9 45#include "AliAODMCParticle.h"
213dd723 46#include "AliAnalysisManager.h"
47#include "AliLog.h"
48#include "TGeoManager.h"
213dd723 49#include "AliMCEventHandler.h"
50#include "AliMCEvent.h"
51#include "AliStack.h"
52#include "AliPHOSGeometry.h"
53#include "AliTriggerAnalysis.h"
54#include "AliEMCALGeometry.h"
92071be9 55#include "AliAnalysisUtils.h"
56#include "AliOADBContainer.h"
213dd723 57
58ClassImp(AliAnalysisTaskTaggedPhotons)
59
60//______________________________________________________________________________
61AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() :
62 AliAnalysisTaskSE(),
63 fPHOSgeom(0x0),
64 fOutputContainer(0x0),
65 fStack(0x0),
66 fTrackEvent(0x0),
67 fPHOSEvent(0x0),
68 fCurrentMixedList(0x0),
69 fTriggerAnalysis(0x0),
92071be9 70 fUtils(0x0),
213dd723 71 fZmax(0.),
72 fZmin(0.),
73 fPhimax(0.),
74 fPhimin(0.),
75 fCentrality(0.),
92071be9 76 fCentBin(0),
b2a1268b 77 fIsMB(0),
78 fIsMC(0)
213dd723 79{
80 //Deafult constructor
81 //no memory allocations
82 for(Int_t i=0;i<10;i++)
83 for(Int_t j=0;j<2;j++)
84 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
92071be9 85 for(Int_t i=0;i<6;i++)
86 fPHOSBadMap[i]=0x0;
213dd723 87}
88//______________________________________________________________________________
89AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) :
90 AliAnalysisTaskSE(name),
91 fPHOSgeom(0x0),
92 fOutputContainer(0x0),
93 fStack(0x0),
94 fTrackEvent(0x0),
95 fPHOSEvent(0x0),
96 fCurrentMixedList(0x0),
97 fTriggerAnalysis(new AliTriggerAnalysis),
92071be9 98 fUtils(0x0),
213dd723 99 fZmax(-60.),
100 fZmin(60.),
101 fPhimax(250.),
102 fPhimin(320.),
103 fCentrality(0.),
92071be9 104 fCentBin(0),
b2a1268b 105 fIsMB(0),
106 fIsMC(0)
213dd723 107{
108 // Constructor.
109
110 // Output slots #0 write into a TH1 container
111 DefineOutput(1,THashList::Class());
112 // Set bad channel map (empty so far)
113 for(Int_t i=0;i<1;i++)
114 for(Int_t j=0;j<5;j++)
115 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
92071be9 116 for(Int_t i=0;i<6;i++)
117 fPHOSBadMap[i]=0x0;
213dd723 118
119
120}
121
122//____________________________________________________________________________
123AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :
124 AliAnalysisTaskSE(ap.GetName()),
125 fPHOSgeom(0x0),
126 fOutputContainer(0x0),
127 fStack(0x0),
128 fTrackEvent(0x0),
129 fPHOSEvent(0x0),
130 fCurrentMixedList(0x0),
131 fTriggerAnalysis(new AliTriggerAnalysis),
92071be9 132 fUtils(0x0),
213dd723 133 fZmax(-60.),
134 fZmin(60.),
135 fPhimax(250.),
136 fPhimin(320.),
137 fCentrality(0.),
92071be9 138 fCentBin(0),
b2a1268b 139 fIsMB(0),
140 fIsMC(0)
213dd723 141{
142 // cpy ctor
143 fZmax=ap.fZmax ;
144 fZmin=ap.fZmin ;
145 fPhimax=ap.fPhimax ;
146 fPhimin=ap.fPhimin ;
147 for(Int_t i=0;i<1;i++)
148 for(Int_t j=0;j<5;j++)
149 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
92071be9 150 for(Int_t i=0;i<6;i++)
151 fPHOSBadMap[i]=0x0;
bffa2e93 152
213dd723 153}
154
155//_____________________________________________________________________________
156AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)
157{
158// assignment operator
159
160 this->~AliAnalysisTaskTaggedPhotons();
161 new(this) AliAnalysisTaskTaggedPhotons(ap);
162 return *this;
163}
164
165//______________________________________________________________________________
166AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()
167{
168 // dtor
169 if(fOutputContainer && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
170 fOutputContainer->Clear() ;
171 delete fOutputContainer ;
172 }
173 for(Int_t i=0;i<10;i++)
174 for(Int_t j=0;j<2;j++)
175 if(fPHOSEvents[i][j]){
176 delete fPHOSEvents[i][j] ;
177 fPHOSEvents[i][j]=0x0 ;
178 }
179}
180//________________________________________________________________________
181void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
182{
183
184
185 // Create the outputs containers
186 if(fOutputContainer != NULL){
187 delete fOutputContainer;
188 }
189 fOutputContainer = new THashList();
190 fOutputContainer->SetOwner(kTRUE);
191 fOutputContainer->SetName(GetName()) ;
192
193 //QA histograms
194 fOutputContainer->Add(new TH1F("hSelEvents","Event selection", 10,0.,10.)) ;
195
196 //vertex distribution
197 fOutputContainer->Add(new TH1F("hNvertexTracks","N of primary tracks from the primary vertex",150,0.,150.));
198 fOutputContainer->Add(new TH1F("hZvertex","Z vertex",200,-50.,+50.));
199 fOutputContainer->Add(new TH2F("hTrackMult","Charged track multiplicity",100,0.,100.,250,0.,500.));
92071be9 200 fOutputContainer->Add(new TH2F("hTrackEtaPhi","Charged track eta vs phi distribution",200,-2.,2.,200,0.,TMath::TwoPi()));
201 fOutputContainer->Add(new TH2F("hTrackEtaPt","Charged track eta vs pt distribution",200,-2.,2.,200,0.,50.));
213dd723 202
203 //centrality
204 fOutputContainer->Add(new TH1F("hCentrality","Ccentrality",100,0.,100.));
42651096 205 fOutputContainer->Add(new TH2F("hPHOSCentrality","PHOS vs centrality",100,0.,100.,100,0.,100.));
213dd723 206 fOutputContainer->Add(new TH2F("hTOF","cluster TOF",200,0.,20.,300,-3.e-6,6.e-6));
92071be9 207 fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,1500.)) ;
213dd723 208
209 fOutputContainer->Add(new TH2F("hCluNXZM1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
210 fOutputContainer->Add(new TH2F("hCluNXZM2","Clu (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
211 fOutputContainer->Add(new TH2F("hCluNXZM3","Clu (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
212 fOutputContainer->Add(new TH2F("hCluEXZM1","Clu E(X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
213 fOutputContainer->Add(new TH2F("hCluEXZM2","Clu E(X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
214 fOutputContainer->Add(new TH2F("hCluEXZM3","Clu E(X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
215
216 fOutputContainer->Add(new TH2F("hCluArea2M1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
217 fOutputContainer->Add(new TH2F("hCluArea2M2","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
218 fOutputContainer->Add(new TH2F("hCluArea2M3","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
219 fOutputContainer->Add(new TH2F("hCluArea3M1","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
220 fOutputContainer->Add(new TH2F("hCluArea3M2","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
221 fOutputContainer->Add(new TH2F("hCluArea3M3","Clu (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
222
92071be9 223 fOutputContainer->Add(new TH2F("hTofM1","TOF in mod1",200,-1.e-6,1.e-6,200,0.,20.)) ;
224 fOutputContainer->Add(new TH2F("hTofM2","TOF in mod2",200,-1.e-6,1.e-6,200,0.,20.)) ;
225 fOutputContainer->Add(new TH2F("hTofM3","TOF in mod3",200,-1.e-6,1.e-6,200,0.,20.)) ;
226
213dd723 227 char cPID[4][5] ;
228 snprintf(cPID[0],5,"All") ;
229 snprintf(cPID[1],5,"Disp");
230 snprintf(cPID[2],5,"CPV") ;
231 snprintf(cPID[3],5,"Both");
232
233
234 const Int_t nPt=400 ;
235 const Double_t ptMax=40. ;
236 const Int_t nM=400 ;
237 const Double_t mMax=1. ;
238
239 const Int_t nCenBin=5;
240 for(Int_t cen=0; cen<nCenBin; cen++){
241
242 for(Int_t iPID=0; iPID<4; iPID++){
243
244 //Inclusive spectra
245 fOutputContainer->Add(new TH1F(Form("hPhot_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
246 fOutputContainer->Add(new TH1F(Form("hPhot_Dist2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
247 fOutputContainer->Add(new TH1F(Form("hPhot_Dist3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
248
249 fOutputContainer->Add(new TH1F(Form("hPhot_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
250 fOutputContainer->Add(new TH1F(Form("hPhot_Area2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
251 fOutputContainer->Add(new TH1F(Form("hPhot_Area3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
252
253 fOutputContainer->Add(new TH1F(Form("hPhot_nStrictTagged_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
254
255 for(Int_t itag=0; itag<18; itag++){
256 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)) ;
257 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)) ;
258 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)) ;
259 }
260 for(Int_t kind=1; kind<33; kind*=2){
261 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)) ;
262 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)) ;
263 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)) ;
264 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)) ;
265 }
266 }
267
268
269 fOutputContainer->Add(new TH1F(Form("hTaggedMult_cent%d",cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
270
271 //Invariant mass distributions for fake corrections
272
273 for(Int_t iPID=0; iPID<4; iPID++){
274 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)) ;
275 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)) ;
276 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)) ;
213dd723 277 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)) ;
278 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)) ;
279 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)) ;
42651096 280
281 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
282 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
283 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
284 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin1_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
285 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin2_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
286 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin3_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
287
288 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
289 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
290 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
291 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin1_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
292 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin2_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
293 fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin3_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
294
295
296
213dd723 297 }
298
299
300 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.)) ;
301 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.)) ;
302 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.)) ;
303 }//centrality
304
305 //MC
306 char partName[4][10] ;
307 snprintf(partName[0],10,"gamma") ;
308 snprintf(partName[1],10,"pi0");
309 snprintf(partName[2],10,"eta") ;
310 snprintf(partName[3],10,"omega");
311
b2a1268b 312 if(fIsMC){
b5f58f9c 313
314 fOutputContainer->Add(new TH1F("hMCConversionRadius","Clusters without label",600,0.,600.)) ;
315 fOutputContainer->Add(new TH2F("hMCRecPi0Vtx","Secondary pi0s",100,0.,10.,600,0.,600.)) ;
316 fOutputContainer->Add(new TH2F("hMCRecEtaVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
317 fOutputContainer->Add(new TH2F("hMCRecOmegaVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
318 fOutputContainer->Add(new TH2F("hMCRecEtaprVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
319 fOutputContainer->Add(new TH2F("hMCRecK0sVtx","Secondary K0s",100,0.,10.,600,0.,600.)) ;
320 fOutputContainer->Add(new TH2F("hMCRecK0lVtx","Secondary K0l",100,0.,10.,600,0.,600.)) ;
321 fOutputContainer->Add(new TH2F("hMCGammaPi0MisConvR","Converted photons",400,0.,40.,600,0.,600.)) ;
322
323 for(Int_t cen=0; cen<nCenBin; cen++){
92071be9 324 for(Int_t ipart=0; ipart<4; ipart++){
325 fOutputContainer->Add(new TH2F(Form("hMCpi0_ptrap_cen%d",cen),"Spectrum of primary photons",100,0.,10.,100,-2.,2.)) ;
326 fOutputContainer->Add(new TH2F(Form("hMCpi0_ptphi_cen%d",cen),"Spectrum of primary photons",100,0.,10.,100,0.,TMath::TwoPi())) ;
327 fOutputContainer->Add(new TH2F(Form("hMC_%s_vertex",partName[ipart]),"vertex",nPt,0.,ptMax,150,0.,600.)) ;
328 fOutputContainer->Add(new TH1F(Form("hMC_all_%s",partName[ipart]),"Spectum (full rapifity)",nPt,0.,ptMax)) ;
329 fOutputContainer->Add(new TH1F(Form("hMC_unitEta_%s",partName[ipart]),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
330 fOutputContainer->Add(new TH1F(Form("hMC_rap_%s",partName[ipart]),"Rapidity",100,-5.,5.)) ;
331 fOutputContainer->Add(new TH1F(Form("hMC_phi_%s",partName[ipart]),"Azimuthal angle",100,0.,TMath::TwoPi())) ;
332 }
213dd723 333
334
92071be9 335 fOutputContainer->Add(new TH2F(Form("LabelsNPrim_cent%d",cen),"vertex",nPt,0.,ptMax,20,0.,20.)) ;
336 fOutputContainer->Add(new TH1F(Form("LabelsGamma_cent%d",cen),"vertex",nPt,0.,ptMax)) ;
337 fOutputContainer->Add(new TH2F(Form("LabelsGammaE_cent%d",cen),"vertex",nPt,0.,ptMax,100,0.,2.)) ;
213dd723 338
92071be9 339 fOutputContainer->Add(new TH1F(Form("hMCRecNoLabel_cent%d",cen),"Clusters without label",nPt,0.,ptMax)) ;
340 fOutputContainer->Add(new TH1F(Form("hMCConversionRadius_cent%d",cen),"Clusters without label",600,0.,600.)) ;
341 fOutputContainer->Add(new TH2F(Form("hMCRecPi0Vtx_cent%d",cen),"Secondary pi0s",100,0.,10.,600,0.,600.)) ;
342 fOutputContainer->Add(new TH2F(Form("hMCRecEtaVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
343 fOutputContainer->Add(new TH2F(Form("hMCRecOmegaVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
344 fOutputContainer->Add(new TH2F(Form("hMCRecEtaprVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
345 fOutputContainer->Add(new TH2F(Form("hMCRecK0sVtx_cent%d",cen),"Secondary K0s",100,0.,10.,600,0.,600.)) ;
346 fOutputContainer->Add(new TH2F(Form("hMCRecK0lVtx_cent%d",cen),"Secondary K0l",100,0.,10.,600,0.,600.)) ;
347 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisConvR_cent%d",cen),"Converted photons",400,0.,40.,600,0.,600.)) ;
213dd723 348
349
92071be9 350 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0PrimMgg_cent%d",cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
351 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0RecMgg_cent%d",cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
213dd723 352
92071be9 353 for(Int_t iPID=0; iPID<4; iPID++){
354 fOutputContainer->Add(new TH1F(Form("hMCRecNoPrim_%s_cent%d",cPID[iPID],cen),"Reconstructed particles withour primary",nPt,0.,ptMax)) ;
355 fOutputContainer->Add(new TH1F(Form("hMCRecGamma_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: Gamma",nPt,0.,ptMax)) ;
356 fOutputContainer->Add(new TH1F(Form("hMCRecPhotConv_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: e+-",nPt,0.,ptMax)) ;
357 fOutputContainer->Add(new TH1F(Form("hMCRecPi0_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: pi0",nPt,0.,ptMax)) ;
358 fOutputContainer->Add(new TH1F(Form("hMCRecEta_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: eta",nPt,0.,ptMax)) ;
359 fOutputContainer->Add(new TH1F(Form("hMCRecOmega_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: Gamma",nPt,0.,ptMax)) ;
360 fOutputContainer->Add(new TH1F(Form("hMCRecEtapr_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: eta prime",nPt,0.,ptMax)) ;
361 fOutputContainer->Add(new TH1F(Form("hMCRecPbar_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: bar(p)",nPt,0.,ptMax)) ;
362 fOutputContainer->Add(new TH1F(Form("hMCRecNbar_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: bar(n)",nPt,0.,ptMax)) ;
363 fOutputContainer->Add(new TH1F(Form("hMCRecPipm_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: pipm",nPt,0.,ptMax)) ;
364 fOutputContainer->Add(new TH1F(Form("hMCRecN_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: n",nPt,0.,ptMax)) ;
365 fOutputContainer->Add(new TH1F(Form("hMCRecP_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: p",nPt,0.,ptMax)) ;
366 fOutputContainer->Add(new TH1F(Form("hMCRecKpm_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K+-",nPt,0.,ptMax)) ;
367 fOutputContainer->Add(new TH1F(Form("hMCRecK0s_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0s",nPt,0.,ptMax)) ;
368 fOutputContainer->Add(new TH1F(Form("hMCRecK0l_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
369 fOutputContainer->Add(new TH1F(Form("hMCRecUnknownCh_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
370 fOutputContainer->Add(new TH1F(Form("hMCRecUnknownNeu_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
371
372 //Decay photons
373 fOutputContainer->Add(new TH1F(Form("hMCRecGammaDir_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, no primary",nPt,0.,ptMax)) ;
374 fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
375 fOutputContainer->Add(new TH1F(Form("hMCRecGammaEta_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from eta",nPt,0.,ptMax)) ;
376 fOutputContainer->Add(new TH1F(Form("hMCRecGammaOmega_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from omega",nPt,0.,ptMax)) ;
377 fOutputContainer->Add(new TH1F(Form("hMCRecGammaOther_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from other",nPt,0.,ptMax)) ;
378 fOutputContainer->Add(new TH1F(Form("hMCRecPhotNoPrim_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax)) ;
213dd723 379
92071be9 380 //Pi0 decay photons
381
382 //MC tagging: reasons of partner loss etc.
383 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_cen%d",cen),"Decay photons with partner not in Stack", nPt,0.,ptMax)) ;
384 fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartn_cen%d",cen),"Decay photon with found partner", nPt,0.,ptMax)) ;
385 fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_cen%d",cen),"Decay photon with wrong mass", nPt,0.,ptMax)) ;
386 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_cen%d",cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax)) ;
387 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_cen%d",cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax)) ;
388 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA2_cen%d",cen),"Decay photons with partner missed due geometry Fid. area. 2", nPt,0.,ptMax)) ;
389 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA3_cen%d",cen),"Decay photons with partner missed due geometry Fid. area. 3", nPt,0.,ptMax)) ;
390 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnConv_cen%d",cen),"Decay photons with partner missed due to conversion", nPt,0.,ptMax)) ;
391 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnEmin_cen%d",cen),"Decay photons with partner missed due to low energy", nPt,0.,ptMax)) ;
392 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnOther_cen%d",cen),"Decay photons with partner missed due unknown reason", nPt,0.,ptMax)) ;
393 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_cen%d",cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax)) ;
394 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_cen%d",cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax)) ;
395
396
397 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_cen%d",cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax)) ;
398 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_cen%d",cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax)) ;
399 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_cen%d",cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax)) ;
400 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_cen%d",cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax)) ;
401 fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_cen%d",cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax)) ;
402 fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_cen%d",cen),"Decay photons with rec partner", nPt,0.,ptMax)) ;
403
404 fOutputContainer->Add(new TH2F(Form("hMCmass_cen%d",cen),"Mass with reconstructed decay partner",nM,0.,mMax,nPt,0.,ptMax)) ;
405
406
407 fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0Dalitz_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
408 fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0NoStack_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
409 fOutputContainer->Add(new TH1F(Form("hMCGammaPi02Gamma_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
410 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0Rec_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
411 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0RecSoft_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
412 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0RecCells_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
b5f58f9c 413
414 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0Tagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
415 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0FakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
416 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0TrueTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
417 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MultyTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
418 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeo_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
419 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA1_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
420 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA2_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
421 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA3_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
422 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisConv_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
423 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisEmin_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
424 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisOther_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
425 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisFakePrim_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
426 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisFoundPrim_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
427 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisNPhot_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
213dd723 428
429 //all clusters fake tagged
b5f58f9c 430 fOutputContainer->Add(new TH1F(Form("hMCAllFakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
213dd723 431
b5f58f9c 432 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)) ;
433 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)) ;
434
213dd723 435 }
b5f58f9c 436 fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisPartner_cent%d",cen),"Spectrum of missed partners",nPt,0.,ptMax)) ;
437 fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisPartnerEtaPhi_cent%d",cen),"Spectrum of missed partners",100,-0.2,0.2,100,4.5,5.6)) ;
438 }
b2a1268b 439 }
440
441 //If we work with MC, need to set Sumw2 - we will use weights
442 if(fIsMC){
443 for(Int_t i=0; i<fOutputContainer->GetSize();i++){
444 ((TH1*)fOutputContainer->At(i))->Sumw2() ;
445 }
446 }
447
448
213dd723 449
450 for(Int_t i=0;i<1;i++)
451 for(Int_t j=0;j<5;j++)
452 fPHOSEvents[i][j]=0x0 ; //Container for PHOS photons
453
454
455 PostData(1, fOutputContainer);
456
457
458}
459
460//______________________________________________________________________________
461void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
462{
463 //Select events
464 //Select photons
465 //Fill QA histograms
466 //Fill Tagging histogsms
92071be9 467
468
469 const Double_t kEcrossCut=0.98 ;
470 const Double_t kTOFMaxCut= 100.e-9 ;
471 const Double_t kTOFMinCut=-100.e-9 ;
213dd723 472
473 // Event selection flags
474
92071be9 475// FillHistogram("hSelEvents",0) ;
213dd723 476
b5f58f9c 477 AliVEvent* event = (AliVEvent*)InputEvent();
213dd723 478 if(!event){
479 AliDebug(1,"No event") ;
480 PostData(1, fOutputContainer);
481 return;
482 }
92071be9 483 FillHistogram("hSelEvents",1) ;
b2a1268b 484
485 //MC stack init
486 fStack = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
213dd723 487
488 //read geometry if not read yet
b2a1268b 489 if(fPHOSgeom==0){
213dd723 490 InitGeometry() ;
b2a1268b 491 }
213dd723 492
92071be9 493 if(!fUtils)
494 fUtils = new AliAnalysisUtils();
495
496 Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7) ;
497 Bool_t isPHI7 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kPHI7);
498
499 if((fIsMB && !isMB) || (fIsMB && isPHI7) || (!fIsMB && !isPHI7)){
500 PostData(1, fOutputContainer);
501 return;
502 }
503 FillHistogram("hSelEvents",2) ;
213dd723 504
505 // Checks if we have a primary vertex
506 // Get primary vertices form AOD
507
213dd723 508 Double_t vtx5[3];
b5f58f9c 509 vtx5[0] = event->GetPrimaryVertex()->GetX();
510 vtx5[1] = event->GetPrimaryVertex()->GetY();
511 vtx5[2] = event->GetPrimaryVertex()->GetZ();
512
513 FillHistogram("hNvertexTracks",event->GetPrimaryVertex()->GetNContributors());
514 FillHistogram("hZvertex" ,vtx5[2]);
515 if (TMath::Abs(vtx5[2]) > 10. ){
213dd723 516 PostData(1, fOutputContainer);
517 return ;
518 }
519
520 FillHistogram("hSelEvents",3) ;
521
522
b5f58f9c 523// if (event->IsPileupFromSPD()){
524// PostData(1, fOutputContainer);
525// return ;
526// }
92071be9 527
528 if(!fUtils->IsVertexSelected2013pA(event)){
529 PostData(1, fOutputContainer);
530 return ;
531 }
213dd723 532 FillHistogram("hSelEvents",4) ;
92071be9 533
534 if(fUtils->IsPileUpEvent(event)){
535 PostData(1, fOutputContainer);
536 return ;
537 }
538 FillHistogram("hSelEvents",5) ;
539
213dd723 540 //centrality
541 AliCentrality *centrality = event->GetCentrality();
542 if( centrality )
543 fCentrality=centrality->GetCentralityPercentile("V0M");
544 else {
545 AliError("Event has 0x0 centrality");
546 fCentrality = -1.;
547 }
548 FillHistogram("hCentrality",fCentrality) ;
549
550 if(fCentrality<0. || fCentrality>=100.){
551 PostData(1, fOutputContainer);
552 return ;
553 }
554 fCentBin = (Int_t)(fCentrality/20.) ;
555
92071be9 556 FillHistogram("hSelEvents",6) ;
213dd723 557
558
559 //Vtx class z-bin
560 Int_t zvtx = 0 ;
561
562 //Calculate charged multiplicity
563 Int_t trackMult = 0;
564 if(fTrackEvent)
565 fTrackEvent->Clear() ;
566 else
0b9a6a23 567 fTrackEvent = new TClonesArray("AliCaloPhoton",event->GetNumberOfTracks()) ;
92071be9 568
213dd723 569 for (Int_t i=0;i<event->GetNumberOfTracks();++i) {
92071be9 570 AliAODTrack *track = (AliAODTrack*)event->GetTrack(i) ;
571 if(!track->IsHybridGlobalConstrainedGlobal())
572 continue ;
213dd723 573 if(TMath::Abs(track->Eta())< 0.9){
574 if(trackMult>=fTrackEvent->GetSize())
575 fTrackEvent->Expand(2*trackMult) ;
0b9a6a23 576 new ((*fTrackEvent)[trackMult]) AliCaloPhoton(track->Px(),track->Py(),track->Pz(),track->P());
213dd723 577 trackMult++;
92071be9 578 FillHistogram("hTrackEtaPhi",track->Eta(),track->Phi()) ;
579 FillHistogram("hTrackEtaPt",track->Eta(),track->Pt()) ;
213dd723 580 }
213dd723 581 }
582 FillHistogram("hTrackMult",fCentrality,trackMult+0.5) ;
583
584 if(!fPHOSEvents[zvtx][fCentBin])
585 fPHOSEvents[zvtx][fCentBin]=new TList() ;
586 fCurrentMixedList = fPHOSEvents[zvtx][fCentBin] ;
587
588 const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
589
590 //---------Select photons-------------------
591 Int_t multClust = event->GetNumberOfCaloClusters();
592 if(!fPHOSEvent)
0b9a6a23 593 fPHOSEvent = new TClonesArray("AliCaloPhoton",multClust);
213dd723 594 else
595 fPHOSEvent->Clear() ;
596 Int_t inList = 0; //counter of caloClusters
597
598 for (Int_t i=0; i<multClust; i++) {
b5f58f9c 599 AliVCluster * clu = event->GetCaloCluster(i);
213dd723 600
601 if(!clu->IsPHOS())
602 continue ;
603
604 if(clu->E()<0.1)
605 continue;
606
607 if(clu->GetNCells()<3)
608 continue ;
609
610 if(clu->GetM02()<0.2)
611 continue ;
42651096 612
92071be9 613 if(clu->GetMCEnergyFraction()>kEcrossCut) //Ecross cut, should be filled with Tender
42651096 614 continue ;
213dd723 615
616 Float_t pos[3] ;
617 clu->GetPosition(pos) ;
618 Int_t fidArea=GetFiducialArea(pos) ;
92071be9 619// if(fidArea==0) //Bad cell
620// continue;
213dd723 621
622 TVector3 global1(pos) ;
623 Int_t relId[4] ;
624 fPHOSgeom->GlobalPos2RelId(global1,relId) ;
625 Int_t mod = relId[0] ;
626 Int_t cellX = relId[2];
627 Int_t cellZ = relId[3] ;
628
92071be9 629 FillHistogram("hTOF",clu->E(),clu->GetTOF()) ;
630 FillHistogram(Form("hTofM%d",mod),clu->GetTOF(),clu->E()) ;
631 if(clu->GetTOF() < kTOFMinCut || clu->GetTOF() > kTOFMaxCut)
632 continue ;
633
634// if(clu->GetDistanceToBadChannel()<2.5)
635// continue ;
636
637
213dd723 638 FillHistogram(Form("hCluNXZM%d",mod),cellX,cellZ,1.);
639 FillHistogram(Form("hCluEXZM%d",mod),cellX,cellZ,clu->E());
640 if(fidArea>1){
641 FillHistogram(Form("hCluArea2M%d",mod),cellX,cellZ,1.);
642 if(fidArea>2){
643 FillHistogram(Form("hCluArea3M%d",mod),cellX,cellZ,1.);
644 }
645 }
646
647 TLorentzVector momentum ;
648 clu->GetMomentum(momentum, vtx5);
0b9a6a23 649 AliCaloPhoton *p = new ((*fPHOSEvent)[inList]) AliCaloPhoton(momentum.Px(),momentum.Py(),momentum.Pz(),clu->E() );
213dd723 650 inList++;
651
652 Int_t isolation = EvalIsolation(&momentum) ;
0b9a6a23 653 p->SetIsolationTag(isolation) ;
213dd723 654
655 p->SetDistToBad((Int_t)(1.+clu->GetDistanceToBadChannel()/2.2));
656
0b9a6a23 657 p->SetTagInfo(0); //Strict PID pi0 partner not found
213dd723 658 p->SetTagged(kFALSE); //Reconstructed pairs found
92071be9 659
213dd723 660 p->SetFiducialArea(fidArea) ;
661
b2a1268b 662 if(fIsMC){
92071be9 663 //This is primary entered PHOS
0b9a6a23 664 FillHistogram(Form("LabelsNPrim_cent%d",fCentBin),clu->E(),float(clu->GetNLabels())) ;
665 Int_t primLabel=clu->GetLabelAt(0) ; //FindPrimary(clu,sure) ;
92071be9 666 //Look what particle left vertex
667 if(primLabel>-1){
668 AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(primLabel) ;
669 Int_t iparent=primLabel;
670 AliAODMCParticle * parent = prim;
671 Double_t r2=prim->Xv()*prim->Xv()+prim->Yv()*prim->Yv() ;
672 while((r2 > rcut*rcut) && (iparent>-1)){
673 iparent=parent->GetMother();
674 parent=(AliAODMCParticle*)fStack->At(iparent);
0b9a6a23 675 r2=parent->Xv()*parent->Xv()+parent->Yv()*parent->Yv() ;
92071be9 676 }
0b9a6a23 677 p->SetPrimary(primLabel) ;
678 p->SetPrimaryAtVertex(iparent) ;
679 p->SetWeight(PrimaryParticleWeight(parent)) ;
92071be9 680 }
681 else{
0b9a6a23 682 p->SetPrimary(-1); //Primary index
683 p->SetPrimaryAtVertex(-1) ;
b2a1268b 684 p->SetWeight(1.) ;
92071be9 685 }
213dd723 686 }
687 else{
0b9a6a23 688 p->SetPrimary(-1); //Primary index
689 p->SetPrimaryAtVertex(-1) ;
b2a1268b 690 p->SetWeight(1.) ;
213dd723 691 }
213dd723 692 //PID criteria
42651096 693 p->SetDispBit(clu->Chi2()<2.5) ;
213dd723 694 p->SetTOFBit(TestTOF(clu->GetTOF(),clu->E())) ;
0b9a6a23 695 p->SetCPVBit(clu->GetEmcCpvDistance()>2.5) ;
213dd723 696 }
42651096 697 FillHistogram("hPHOSCentrality",fCentrality,inList+0.5) ;
213dd723 698
699 FillMCHistos() ;
700 FillTaggingHistos() ;
701
702 //Remove old events
703 fCurrentMixedList->AddFirst(fPHOSEvent);
704 fPHOSEvent=0x0 ;
705 if(fCurrentMixedList->GetSize() > 10){
706 TClonesArray *tmp = static_cast <TClonesArray*> (fCurrentMixedList->Last());
707 fCurrentMixedList->Remove(tmp);
708 delete tmp;
709 }
710
711 PostData(1, fOutputContainer);
712
713}
714//________________________________________________
715void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
716
717 //MC info about this particle
b2a1268b 718 if(!fIsMC)
213dd723 719 return ;
213dd723 720 const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
92071be9 721 const Double_t phiMin=260.*TMath::Pi()/180. ;
722 const Double_t phiMax=320.*TMath::Pi()/180. ;
723
724 Int_t nPrim = fStack->GetEntriesFast() ;
725 //Fill Primary particl yields
213dd723 726
92071be9 727 for(Int_t i=0;i<nPrim;i++){
728 AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(i) ;
729 Double_t r2=prim->Xv()*prim->Xv()+prim->Yv()*prim->Yv() ;
730 if(r2>rcut*rcut)
731 continue ;
732
733 Int_t pdg=prim->GetPdgCode() ;
734 char partName[30] ;
735 if(pdg == 111)
736 snprintf(partName,30,"pi0") ;
737 else
738 if(pdg == 221)
739 snprintf(partName,30,"eta") ;
740 else
741 if(pdg == 22)
742 snprintf(partName,30,"gamma") ;
743 else
744 if(pdg == 310)
745 snprintf(partName,30,"K0s") ;
746 else
747 if(abs(pdg) == 321)
748 snprintf(partName,30,"Kpm") ;
749 else
750 if(abs(pdg) == 211)
751 snprintf(partName,30,"pipm") ;
752 else
753 continue ;
754
755
213dd723 756 //Primary particle
92071be9 757 Double_t phi=prim->Phi() ;
758 while(phi<0.)phi+=TMath::TwoPi() ;
759 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
760 Double_t pt=prim->Pt() ;
213dd723 761
213dd723 762 //Total number of pi0 with creation radius <1 cm
92071be9 763 Double_t w = PrimaryParticleWeight(prim) ;
213dd723 764 FillHistogram(Form("hMC_all_%s",partName),pt,w) ;
92071be9 765 if(TMath::Abs(prim->Y())<0.13){
766 FillHistogram(Form("hMC_phi_%s",partName),phi,w) ;
767 if(phi>phiMin && phi<phiMax)
768 FillHistogram(Form("hMC_%s_pt_cen%d",partName,fCentBin),pt,w) ;
213dd723 769 }
770
92071be9 771 FillHistogram(Form("hMC_rap_%s",partName),prim->Y(),w) ;
772 //Some additional QA
773 if(pdg == 111){
774 FillHistogram(Form("hMCpi0_ptrap_cen%d",fCentBin),pt,prim->Y(),w) ;
775 FillHistogram(Form("hMCpi0_ptphi_cen%d",fCentBin),pt,phi,w) ;
776 }
777 if(pdg == 22){
778 FillHistogram(Form("hMCgamma_ptrap_cen%d",fCentBin),pt,prim->Y(),w) ;
779 FillHistogram(Form("hMCgamma_ptphi_cen%d",fCentBin),pt,phi,w) ;
780 }
213dd723 781
213dd723 782 }
92071be9 783
213dd723 784
785
786 //Clussify reconstructed clusters
787 //First - photons (from vertex) and contaminations
788 //Second - photons from different sources
789 //Third - photons from pi0s - missed for different reasons
790
791 const Int_t n=fPHOSEvent->GetEntriesFast() ;
792 for(Int_t i=0;i<n;i++){
0b9a6a23 793 AliCaloPhoton *p = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
794 Int_t label=p->GetPrimary() ;
92071be9 795 if(label<0){ //No label!
b2a1268b 796 FillHistogram("hMCRecNoLabel",p->Pt(),p->GetWeight());
213dd723 797 continue ;
798 }
213dd723 799
92071be9 800
0b9a6a23 801 AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(p->GetPrimary()) ;
92071be9 802 //Look what particle left virtex
0b9a6a23 803 Int_t iparent=p->GetPrimary();
92071be9 804 AliAODMCParticle * parent = prim;
805 while(parent->Xv()*parent->Xv()+parent->Yv()*parent->Yv() > rcut*rcut){
806 iparent=parent->GetMother();
807 if(iparent<0)
213dd723 808 break ;
92071be9 809 parent = (AliAODMCParticle*)fStack->At(iparent) ;
810 }
811 Int_t parentPDG=parent->GetPdgCode() ;
812 switch(parentPDG){
813 case 22: //electron/positron conversion
814 FillPIDHistograms("hMCRecPhoton",p); //Reconstructed with photon from conversion primary
213dd723 815 break ;
92071be9 816 case 11:
817 case -11: //electron/positron conversion
818 FillPIDHistograms("hMCRecE",p); //Reconstructed with photon from conversion primary
213dd723 819 break ;
92071be9 820 case -2212:
213dd723 821 FillPIDHistograms("hMCRecPbar",p); //Reconstructed with photon from antibaryon annihilation
822 break ;
92071be9 823 case -2112: //antineutron & antiproton conversion
213dd723 824 FillPIDHistograms("hMCRecNbar",p); //Reconstructed with photon from antibaryon annihilation
825 break ;
92071be9 826 case 211:
827 case -211:
828 case 2212:
829 case 321:
830 case -321:
831 FillPIDHistograms("hMCRecCharg",p); //Reconstructed with photon from conversion primary
213dd723 832 break ;
92071be9 833 case 310:
213dd723 834 FillPIDHistograms("hMCRecK0s",p); //Reconstructed with photon from conversion primary
835 break ;
92071be9 836 case 2112: //antineutron & antiproton conversion
837 case 130:
838 FillPIDHistograms("hMCRecNeutral",p); //Reconstructed with photon from antibaryon annihilation
213dd723 839 break ;
92071be9 840 case -1: //direct photon or no primary
841 FillPIDHistograms("hMCRecNoPRim",p);
842 break ;
843 default:
844 printf("Unknown PDG: %d \n",parentPDG) ;
845 FillPIDHistograms("hMCRecUnknown",p);
846 break ;
847 }
213dd723 848
92071be9 849
850 //Now classify decay photon
851 if(parentPDG==22){
852 Int_t iGrandParent=parent->GetMother();
853 if(iGrandParent<0 || iGrandParent>=fStack->GetEntriesFast()){
854 FillPIDHistograms("hMCRecPhotNoPrim",p);
855 continue ;
213dd723 856 }
92071be9 857 AliAODMCParticle * grandParent = (AliAODMCParticle*)fStack->At(iGrandParent) ;
858 Int_t grandParentPDG=grandParent->GetPdgCode() ;
859 switch(grandParentPDG){
860 case 111: //pi0
861 FillPIDHistograms("hMCRecPhotPi0",p);
862 break ;
863 case 221: //eta decay
864 FillPIDHistograms("hMCRecPhotEta",p);
865 break ;
866 case 223: //omega meson decay
867 FillPIDHistograms("hMCRecPhotOmega",p);
868 break ;
869 default:
870 FillPIDHistograms("hMCRecPhotOther",p);
213dd723 871 }
92071be9 872 //--------consider pi0 decays--------------------
873 if(grandParentPDG==111){
874 //First find which daughter is our cluster
875 //iparent - index of curent photon
876 Int_t ipartner=grandParent->GetDaughter(0) ;
877 if(ipartner==iparent){//look for other
878 if(grandParent->GetNDaughters()>1){
879 ipartner=grandParent->GetDaughter(1);
213dd723 880 }
881 else{
92071be9 882 ipartner=-1 ;
213dd723 883 }
884 }
92071be9 885
886 //There is no partner in stack
887 if(ipartner==-1){
888 FillPIDHistograms("hMCDecWMisPartnStack",p) ;
889 }
890 else{
891 AliAODMCParticle * partner = (AliAODMCParticle *)fStack->At(ipartner);
892 //Check if partner is registered and made correct mass
893 //If not - trace the reason
0b9a6a23 894 AliCaloPhoton *pp = 0x0 ;
213dd723 895
92071be9 896 for(Int_t ii=0;ii<n;ii++){
897 if(i==ii) continue;
0b9a6a23 898 AliCaloPhoton * tmp = static_cast<AliCaloPhoton*>(fPHOSEvent->At(ii));
899 Int_t ipartnPrim = tmp->GetPrimary() ;
92071be9 900 while(ipartnPrim>-1){
901 if(ipartnPrim==ipartner)
902 break ;
903 ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
904 }
905 if(ipartnPrim==ipartner){
906 pp=tmp ;
907 break ;
908 }
909 }
910
911 if(pp){
912 //Partner reconstructed, but did not pass cuts
92071be9 913 FillPIDHistograms("hMCDecWRecPartn",p) ;
0b9a6a23 914 Double_t invMass=(*p+ *pp).M() ;
b2a1268b 915 FillHistogram("hMCmass",invMass,p->Pt(),p->GetWeight()) ;
92071be9 916 if(IsInPi0Band(invMass,p->Pt())){
917 FillPIDHistograms("hMCDecWithFoundPartn",p) ;
918 }
919 else{
920 FillPIDHistograms("hMCDecWithWrongMass",p) ;
921 }
92071be9 922 }
923 else{//Partner not reconstructed
924 if(partner->GetPdgCode()==22){
213dd723 925 Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
926
92071be9 927 //Check if partner miss PHOS
928 Int_t modulenum ;
213dd723 929 Double_t ztmp=0.,xtmp=0. ;
92071be9 930 Double_t vtx[3]={partner->Xv(),partner->Yv(),partner->Zv()} ;
931 Bool_t impact=fPHOSgeom->ImpactOnEmc(vtx,partner->Theta(),partner->Phi(),modulenum,ztmp,xtmp) ;
932
933 if(impact){//still check bad map
934 Int_t relId[4] ;
935 fPHOSgeom->RelPosToRelId(modulenum,xtmp,ztmp,relId) ;
936 if ( !IsGoodChannel(modulenum,relId[2],relId[3]) ) {
937 impact=kFALSE ;
938 }
939 }
940
941 if(!impact){ //this photon cannot hit PHOS
942 FillPIDHistograms("hMCDecWMisPartnAccept",p) ; //Spectrum of tagged with missed partner
213dd723 943 Int_t iFidArea = p->GetFiducialArea();
944 if(iFidArea>0){
92071be9 945 FillPIDHistograms("hMCDecWMisPartnAcceptFA1",p) ; //Spectrum of tagged with missed partner
213dd723 946 if(iFidArea>1){
92071be9 947 FillPIDHistograms("hMCDecWMisPartnAcceptFA2",p) ; //Spectrum of tagged with missed partner
213dd723 948 if(iFidArea>2){
92071be9 949 FillPIDHistograms("hMCDecWMisPartnAcceptFA3",p) ; //Spectrum of tagged with missed partner
213dd723 950 }
951 }
952 }
953 isPartnerLost=kTRUE;
954 }
213dd723 955
92071be9 956 if(!isPartnerLost){
957 //this photon is converted before it is registered
958 if(partner->GetNDaughters()>0){
959 AliAODMCParticle* tmpP=(AliAODMCParticle*)fStack->At(partner->GetDaughter(0));
960 if(tmpP->Xv()*tmpP->Xv()+tmpP->Yv()*tmpP->Yv()<450.*450.){
961 FillPIDHistograms("hMCDecWMisPartnConv",p) ; //Spectrum of tagged with missed partner
962 isPartnerLost=kTRUE;
963 }
964 }
965 }
213dd723 966 if(!isPartnerLost &&
92071be9 967 partner->E()<0.3){ //energy is not enough to be registered by PHOS
968 FillPIDHistograms("hMCDecWMisPartnEmin",p) ; //Spectrum of tagged with missed partner
213dd723 969 isPartnerLost=kTRUE;
970 }
92071be9 971 if(!isPartnerLost){ //Reason not found!!!!!
972 FillPIDHistograms("hMCDecWMisPartnOther",p);
973 }
974 else{//Sum of all missed partners
975 FillPIDHistograms("hMCDecWMisPartnAll",p);
213dd723 976 }
977 }//Partner - photon
978 else{//partner not photon
92071be9 979 FillPIDHistograms("hMCDecWMisPartnNPhot",p);
213dd723 980 }
92071be9 981
982 }//Partner not reconstructed
983 }//Partner in stack
984 }//photon from pi0 decay
985 }//photon
986 } //PHOS clusters
213dd723 987
988}
92071be9 989
213dd723 990//________________________________________________
991void AliAnalysisTaskTaggedPhotons::FillTaggingHistos(){
992 //Fill all necessary histograms
993
994
995 //First fill invariant mass distrubutions and mark tagged photons
996 //Invariant Mass analysis
997 const Int_t n=fPHOSEvent->GetEntriesFast() ;
998 for(Int_t i=0;i<n-1;i++){
0b9a6a23 999 AliCaloPhoton *p1 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
213dd723 1000 for(Int_t j = i+1 ; j < n ; j++) {
0b9a6a23 1001 AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(j));
213dd723 1002
0b9a6a23 1003 Double_t invMass = (*p1 + *p2).M();
42651096 1004
1005 if((p1->E()>0.1) && (p2->E()>0.1)){
1006 FillPIDHistograms("hInvM_Re_Emin1",p1,p2,invMass) ;
1007 if((p1->E()>0.2) && (p2->E()>0.2)){
1008 FillPIDHistograms("hInvM_Re_Emin2",p1,p2,invMass) ;
1009 if((p1->E()>0.3) && (p2->E()>0.3)){
1010 FillPIDHistograms("hInvM_Re_Emin3",p1,p2,invMass) ;
1011 }
1012 }
1013 }
1014
213dd723 1015
1016 if(p2->E()>0.1){
42651096 1017 FillPIDHistograms("hSingleInvM_Re_Emin1",p1,invMass) ;
213dd723 1018 if(p2->E()>0.2){
42651096 1019 FillPIDHistograms("hSingleInvM_Re_Emin2",p1,invMass) ;
213dd723 1020 if(p2->E()>0.3){
42651096 1021 FillPIDHistograms("hSingleInvM_Re_Emin3",p1,invMass) ;
213dd723 1022 }
1023 }
1024 }
1025
1026 if(p1->E()>0.1){
42651096 1027 FillPIDHistograms("hSingleInvM_Re_Emin1",p2,invMass) ;
213dd723 1028 if(p1->E()>0.2){
42651096 1029 FillPIDHistograms("hSingleInvM_Re_Emin2",p2,invMass) ;
213dd723 1030 if(p1->E()>0.3){
42651096 1031 FillPIDHistograms("hSingleInvM_Re_Emin3",p2,invMass) ;
213dd723 1032 }
1033 }
1034 }
1035 if(TestPID(3, p2)){
1036 if(p2->E()>0.1){
42651096 1037 FillPIDHistograms("hSingleInvM_Re_Emin1_Strict",p1,invMass) ;
213dd723 1038 if(p2->E()>0.2){
42651096 1039 FillPIDHistograms("hSingleInvM_Re_Emin2_Strict",p1,invMass) ;
213dd723 1040 if(p2->E()>0.3){
42651096 1041 FillPIDHistograms("hSingleInvM_Re_Emin3_Strict",p1,invMass) ;
213dd723 1042 }
1043 }
1044 }
1045 }
1046 if(TestPID(3, p1)){
1047 if(p1->E()>0.1){
42651096 1048 FillPIDHistograms("hSingleInvM_Re_Emin1_Strict",p2,invMass) ;
213dd723 1049 if(p1->E()>0.1){
42651096 1050 FillPIDHistograms("hSingleInvM_Re_Emin2_Strict",p2,invMass) ;
213dd723 1051 if(p1->E()>0.3){
42651096 1052 FillPIDHistograms("hSingleInvM_Re_Emin3_Strict",p2,invMass) ;
213dd723 1053 }
1054 }
1055 }
1056 }
b2a1268b 1057 if(IsSameParent(p1,p2)==111){
213dd723 1058 FillPIDHistograms("hMC_InvM_Re",p1,invMass) ;
1059 FillPIDHistograms("hMC_InvM_Re",p2,invMass) ;
1060 if(TestPID(3, p2)){
1061 FillPIDHistograms("hMC_InvM_Re_Strict",p1,invMass) ;
1062 }
1063 if(TestPID(3, p1)){
1064 FillPIDHistograms("hMC_InvM_Re_Strict",p2,invMass) ;
1065 }
1066 }
1067
1068 //Tagging: 1,2,3 sigma
1069 //Emin=100,200,300 Mev
1070 //IsInPi0Band(mass, Ptphoton, type Emin cut
1071 Int_t tag1=0 ;
1072 for(Int_t eminType=0; eminType<3; eminType++){
1073 if(p2->E()>0.1*(eminType+1)){
1074 //Set corresponding bit
92071be9 1075 Double_t nsigma = IsInPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
213dd723 1076 for(Int_t isigma=0; isigma<3; isigma++){
1077 if(nsigma<1+isigma){
1078 tag1|= (1 << (3*eminType+isigma)) ;
1079 if(TestPID(3, p2))
1080 tag1|= (1 << (3*eminType+isigma+9)) ;
1081 }
1082 }
1083 }
1084 }
0b9a6a23 1085 p1->SetTagInfo(tag1) ;
213dd723 1086 Int_t tag2=0 ;
1087 for(Int_t eminType=0; eminType<3; eminType++){
1088 if(p1->E()>0.1*(eminType+1)){
1089 //Set corresponding bit
92071be9 1090 Double_t nsigma = IsInPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
213dd723 1091 for(Int_t isigma=0; isigma<3; isigma++){
1092 if(nsigma<1+isigma){
1093 tag2|= (1 << (3*eminType+isigma)) ;
1094 if(TestPID(3, p2))
1095 tag2|= (1 << (3*eminType+isigma+9)) ;
1096 }
1097 }
1098 }
1099 }
0b9a6a23 1100 p2->SetTagInfo(tag2) ;
213dd723 1101
1102 if(tag1 & (1<<7)){ //2 sigma, Emin=0.3: default tagging
1103 if(p1->IsTagged()){//Multiple tagging
b2a1268b 1104 FillHistogram(Form("hTaggedMult_cent%d",fCentBin),p1->Pt(),p1->GetWeight());
213dd723 1105 }
1106 p1->SetTagged(kTRUE) ;
1107 }
1108 if(tag2 & (1<<7)){ //2 sigma, Emin=0.3: default tagging
1109 if(p2->IsTagged()){//Multiple tagging
b2a1268b 1110 FillHistogram(Form("hTaggedMult_cent%d",fCentBin),p2->Pt(),p2->GetWeight());
213dd723 1111 }
1112 p2->SetTagged(kTRUE) ;
1113 }
1114 }
1115 }
1116
1117
1118
1119 //Single particle histgams
1120 for(Int_t i=0;i<n;i++){
0b9a6a23 1121 AliCaloPhoton *p = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
213dd723 1122
0b9a6a23 1123 Int_t isolation = p->GetIsolationTag();
213dd723 1124
1125 //Inclusive spectra
1126 FillPIDHistograms("hPhot",p) ;
1127
1128 if(p->DistToBad()>1){
1129 FillPIDHistograms("hPhot_Dist2",p) ;
1130 if(p->DistToBad()>2){
1131 FillPIDHistograms("hPhot_Dist3",p) ;
1132 }
1133 }
1134
1135
1136 for(Int_t kind=1; kind<33; kind*=2){
1137 if((isolation&kind)){
1138 FillPIDHistograms(Form("hPhot_Isolation%d",kind),p) ;
1139 }
1140 }
1141
1142 Int_t iFidArea = p->GetFiducialArea();
1143 if(iFidArea>0){
1144 FillPIDHistograms("hPhot_Area1",p) ;
1145 for(Int_t kind=1; kind<33; kind*=2){
1146 if((isolation&kind)){
1147 FillPIDHistograms(Form("hPhot_Isolation%d_Area1",kind),p) ;
1148 }
1149 }
1150
1151 //Fill taggings with
1152 //3 Emin cuts
1153 //Default Emin, 1,2,3 sigmas
1154 //strict and loose PID cut on partner
0b9a6a23 1155 Int_t tag=p->GetTagInfo() ;
213dd723 1156 for(Int_t ibit=0; ibit<18; ibit++){
1157 if((tag & (1<<ibit))==0){
1158 FillPIDHistograms(Form("hPhot_nTagged%d_Area1",ibit),p) ;
1159// for(Int_t kind=1; kind<33; kind*=2){
1160// if((isolation&kind)){
1161// FillPIDHistograms(Form("hPhot_nStrictTagged_Isolation%d_Area1",kind),p) ;
1162// }
1163// }
1164 }
1165 }
1166
1167 if(iFidArea>1){
1168 FillPIDHistograms("hPhot_Area2",p) ;
1169 for(Int_t ibit=0; ibit<18; ibit++){
1170 if((tag & (1<<ibit))==0){
1171 FillPIDHistograms(Form("hPhot_nTagged%d_Area2",ibit),p) ;
1172 }
1173 }
1174 if(iFidArea>2){
1175 FillPIDHistograms("hPhot_Area3",p) ;
1176 for(Int_t ibit=0; ibit<18; ibit++){
1177 if((tag & (1<<ibit))==0){
1178 FillPIDHistograms(Form("hPhot_nTagged%d_Area3",ibit),p) ;
1179 }
1180 }
1181 }
1182 }
1183 }
1184 }
1185
1186 //Fill Mixed InvMass distributions:
1187 TIter nextEv(fCurrentMixedList) ;
1188 for(Int_t i=0;i<n;i++){
0b9a6a23 1189 AliCaloPhoton *p1 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
213dd723 1190 while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
1191 Int_t nPhotons2 = event2->GetEntriesFast() ;
1192 for(Int_t j=0; j < nPhotons2 ; j++){
0b9a6a23 1193 AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(event2->At(j)) ;
1194 Double_t invMass = (*p1 + *p2).M();
213dd723 1195
42651096 1196 if((p1->E()>0.1) && (p2->E()>0.1)){
1197 FillPIDHistograms("hInvM_Mi_Emin1",p1,p2,invMass) ;
1198 if((p1->E())>0.2 && (p2->E()>0.2)){
1199 FillPIDHistograms("hInvM_Mi_Emin2",p1,p2,invMass) ;
1200 if((p1->E())>0.3 && (p2->E()>0.3)){
1201 FillPIDHistograms("hInvM_Mi_Emin3",p1,p2,invMass) ;
1202 }
1203 }
1204 }
1205
1206
213dd723 1207 if(p2->E()>0.1){
42651096 1208 FillPIDHistograms("hSingleInvM_Mi_Emin1",p1,invMass) ;
213dd723 1209 if(p2->E()>0.2){
42651096 1210 FillPIDHistograms("hSingleInvM_Mi_Emin2",p1,invMass) ;
213dd723 1211 if(p2->E()>0.3){
42651096 1212 FillPIDHistograms("hSingleInvM_Mi_Emin3",p1,invMass) ;
213dd723 1213 }
1214 }
1215 }
1216 if(TestPID(3, p2)){
1217 if(p2->E()>0.1){
42651096 1218 FillPIDHistograms("hSingleInvM_Mi_Emin1_Strict",p1,invMass) ;
213dd723 1219 if(p2->E()>0.2){
42651096 1220 FillPIDHistograms("hSingleInvM_Mi_Emin2_Strict",p1,invMass) ;
213dd723 1221 if(p2->E()>0.3){
42651096 1222 FillPIDHistograms("hSingleInvM_Mi_Emin3_Strict",p1,invMass) ;
213dd723 1223 }
1224 }
1225 }
1226 }
1227
1228 if(p1->E()>0.1){
42651096 1229 FillPIDHistograms("hSingleInvM_Mi_Emin1",p2,invMass) ;
213dd723 1230 if(p1->E()>0.2){
42651096 1231 FillPIDHistograms("hSingleInvM_Mi_Emin2",p2,invMass) ;
213dd723 1232 if(p1->E()>0.3){
42651096 1233 FillPIDHistograms("hSingleInvM_Mi_Emin3",p2,invMass) ;
213dd723 1234 }
1235 }
1236 }
1237 if(TestPID(3, p1)){
1238 if(p1->E()>0.1){
42651096 1239 FillPIDHistograms("hSingleInvM_Mi_Emin1_Strict",p2,invMass) ;
213dd723 1240 if(p1->E()>0.2){
42651096 1241 FillPIDHistograms("hSingleInvM_Mi_Emin2_Strict",p2,invMass) ;
213dd723 1242 if(p1->E()>0.3){
42651096 1243 FillPIDHistograms("hSingleInvM_Mi_Emin3_Strict",p2,invMass) ;
213dd723 1244 }
1245 }
1246 }
1247 }
1248 }
1249 }
1250 }
213dd723 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//______________________________________________________________________________
92071be9 1267Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
213dd723 1268{
1269 //Parameterization of the pi0 peak region
42651096 1270 //for LHC13bcdef
1271 Double_t mpi0mean = 0.13447 - 1.41259e-03 * TMath::Exp(-pt/1.30044) ;
213dd723 1272
92071be9 1273 Double_t mpi0sigma=TMath::Sqrt(5.22245e-03*5.22245e-03 +2.86851e-03*2.86851e-03/pt) + 9.09932e-05*pt ;
213dd723 1274
1275 return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
1276}
1277//______________________________________________________________________________
b2a1268b 1278Int_t AliAnalysisTaskTaggedPhotons::IsSameParent(const AliCaloPhoton *p1, const AliCaloPhoton *p2)const{
213dd723 1279 //Looks through parents and finds if there was commont pi0 among ancestors
1280
b2a1268b 1281 if(!fIsMC)
1282 return 0 ; //can not say anything
213dd723 1283
0b9a6a23 1284 Int_t prim1 = p1->GetPrimary();
213dd723 1285 while(prim1!=-1){
0b9a6a23 1286 Int_t prim2 = p2->GetPrimary();
213dd723 1287 while(prim2!=-1){
1288 if(prim1==prim2){
b2a1268b 1289 return ((AliAODMCParticle*)fStack->At(prim1))->GetPdgCode() ;
213dd723 1290 }
92071be9 1291 prim2=((AliAODMCParticle*)fStack->At(prim2))->GetMother() ;
213dd723 1292 }
92071be9 1293 prim1=((AliAODMCParticle*)fStack->At(prim1))->GetMother() ;
213dd723 1294 }
b2a1268b 1295 return 0 ;
213dd723 1296}
1297//______________________________________________________________________________
1298Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(const Float_t * position)const{
1299 //calculates in which kind of fiducial area photon hit
1300
1301 TVector3 global1(position) ;
1302 Int_t relId[4] ;
1303 fPHOSgeom->GlobalPos2RelId(global1,relId) ;
92071be9 1304// Int_t mod = relId[0] ;
213dd723 1305 Int_t cellX = relId[2];
1306 Int_t cellZ = relId[3] ;
1307
213dd723 1308 //New we are in good channel,
1309 //calculate distance to the closest group of bad channels
1310 const Int_t cut1=10;
1311 const Int_t cut2=15;
1312//For 3-mod configuration
1313// if((mod==3 && cellX<=cut1) || (mod==1 && cellX>=65-cut1) || cellZ<=cut1 || cellZ>=57-cut1)
1314//For 1+3 configuraion
1315 if( cellX<=cut1 || cellX>=65-cut1 || cellZ<=cut1 || cellZ>=57-cut1)
1316 return 1;
1317// //and from large dead area
213dd723 1318//Full configuration
1319// if((mod==3 && cellX<=cut2) || (mod==1 && cellX>=65-cut2) || cellZ<=cut2 || cellZ>=57-cut2)
1320//1+3 configuration
1321 if( cellX<=cut2 || cellX>=65-cut2 || cellZ<=cut2 || cellZ>=57-cut2)
1322 return 2;
213dd723 1323 //Very good channel
1324 return 3 ;
1325
1326}
1327//______________________________________________________________________________^M
1328void AliAnalysisTaskTaggedPhotons::InitGeometry(){
1329 //Rotation matrixes are set with Tender
1330
1331 if(fPHOSgeom) return ;
213dd723 1332
92071be9 1333
1334 fPHOSgeom = AliPHOSGeometry::GetInstance() ;
1335
1336 if(!fPHOSgeom){ //Geometry not yet constructed with Tender
1337 fPHOSgeom = AliPHOSGeometry::GetInstance("IHEP","");
1338
1339 AliOADBContainer geomContainer("phosGeo");
1340 geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
1341 TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(170000,"PHOSRotationMatrixes");
1342 for(Int_t mod=0; mod<5; mod++) {
1343 if(!matrixes->At(mod)) continue;
1344 fPHOSgeom->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
1345 }
1346 }
1347
1348 //Read BadMap for MC simulations
1349 Int_t runNumber=196208 ; //LHC13bcdef
1350 AliOADBContainer badmapContainer(Form("phosBadMap"));
1351 badmapContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root","phosBadMap");
1352 TObjArray *maps = (TObjArray*)badmapContainer.GetObject(runNumber,"phosBadMap");
1353 if(!maps){
1354 AliError("TaggedPhotons: Can not read Bad map\n") ;
1355 }
1356 else{
1357 AliInfo(Form("TaggedPhotons: Setting PHOS bad map with name %s \n",maps->GetName())) ;
1358 for(Int_t mod=0; mod<5;mod++){
1359 if(fPHOSBadMap[mod])
1360 delete fPHOSBadMap[mod] ;
1361 TH2I * h = (TH2I*)maps->At(mod) ;
1362 if(h)
1363 fPHOSBadMap[mod]=new TH2I(*h) ;
1364 }
1365 }
213dd723 1366}
1367//_____________________________________________________________________________
1368void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x)const{
1369 //FillHistogram
1370 TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
1371 if(tmpI){
1372 tmpI->Fill(x) ;
1373 return ;
1374 }
1375 TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
1376 if(tmpF){
1377 tmpF->Fill(x) ;
1378 return ;
1379 }
1380 TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
1381 if(tmpD){
1382 tmpD->Fill(x) ;
1383 return ;
1384 }
1385 AliInfo(Form("can not find histogram <%s> ",key)) ;
1386}
1387//_____________________________________________________________________________
1388void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x,Double_t y)const{
1389 //FillHistogram
1390 TObject * tmp = fOutputContainer->FindObject(key) ;
1391 if(!tmp){
1392 AliInfo(Form("can not find histogram <%s> ",key)) ;
1393 return ;
1394 }
1395 if(tmp->IsA() == TClass::GetClass("TH1F")){
1396 ((TH1F*)tmp)->Fill(x,y) ;
1397 return ;
1398 }
1399 if(tmp->IsA() == TClass::GetClass("TH2F")){
1400 ((TH2F*)tmp)->Fill(x,y) ;
1401 return ;
1402 }
1403 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1404}
1405
1406//_____________________________________________________________________________
1407void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1408 //Fills 1D histograms with key
1409 TObject * tmp = fOutputContainer->FindObject(key) ;
1410 if(!tmp){
1411 AliInfo(Form("can not find histogram <%s> ",key)) ;
1412 return ;
1413 }
1414 if(tmp->IsA() == TClass::GetClass("TH2F")){
1415 ((TH2F*)tmp)->Fill(x,y,z) ;
1416 return ;
1417 }
1418 if(tmp->IsA() == TClass::GetClass("TH3F")){
1419 ((TH3F*)tmp)->Fill(x,y,z) ;
1420 return ;
1421 }
1422}
1423//_____________________________________________________________________________
0b9a6a23 1424void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p) const{
213dd723 1425
b2a1268b 1426 FillHistogram(Form("%s_All_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
0b9a6a23 1427 if(p->IsDispOK())
b2a1268b 1428 FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
0b9a6a23 1429 if(p->IsCPVOK())
b2a1268b 1430 FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
0b9a6a23 1431 if(p->IsDispOK() && p->IsCPVOK())
b2a1268b 1432 FillHistogram(Form("%s_Both_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
213dd723 1433
1434}
1435//_____________________________________________________________________________
0b9a6a23 1436void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p,Double_t x) const{
213dd723 1437
b2a1268b 1438 FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
0b9a6a23 1439 if(p->IsDispOK())
b2a1268b 1440 FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
0b9a6a23 1441 if(p->IsCPVOK())
b2a1268b 1442 FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
0b9a6a23 1443 if(p->IsDispOK() && p->IsCPVOK())
b2a1268b 1444 FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
213dd723 1445
42651096 1446}
1447//_____________________________________________________________________________
0b9a6a23 1448void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p1,const AliCaloPhoton * p2,Double_t x) const{
42651096 1449
0b9a6a23 1450 Double_t ptPi = (*p1 + *p2).Pt() ;
b2a1268b 1451 Double_t w=p1->GetWeight()*p2->GetWeight() ;
1452 if(IsSameParent(p1,p2)>0)
1453 w=p1->GetWeight() ;
1454 FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,ptPi,w) ;
0b9a6a23 1455 if(p1->IsDispOK() && p2->IsDispOK())
b2a1268b 1456 FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,ptPi,w) ;
0b9a6a23 1457 if(p1->IsCPVOK() && p2->IsCPVOK())
b2a1268b 1458 FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,ptPi,w) ;
0b9a6a23 1459 if(p1->IsDispOK() && p1->IsCPVOK() && p2->IsDispOK() && p2->IsCPVOK())
b2a1268b 1460 FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,ptPi,w) ;
42651096 1461
213dd723 1462}
1463//_____________________________________________________________________________
1464Bool_t AliAnalysisTaskTaggedPhotons::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1465
1466 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1467 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1468 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1469 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1470 Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1471 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1472 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1473 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1474 return (R2<2.5*2.5) ;
1475
1476}
1477//_________________________________________________________________________________
1478Int_t AliAnalysisTaskTaggedPhotons::EvalIsolation(TLorentzVector * ph){
1479
1480 // Check if this particle is isolated.
1481 //We use several cone radii and epsilons.
1482 //As well we look at charged particles and EMCAL ones
1483
1484 const Double_t coneR1=0.3 ;
1485 const Double_t coneR2=0.4 ;
1486 const Double_t coneR3=0.5 ;
1487
1488 const Double_t epsilon1=0.1 ;
1489 const Double_t epsilon2=0.05 ;
1490
1491 if(!ph) return 0 ;
1492
1493 //Sum of energies in cones, tracks and clusters in EMCAL
1494 Double_t eCone1 = 0;
1495 Double_t eCone2 = 0;
1496 Double_t eCone3 = 0;
1497
1498 Double_t phiTrig = ph->Phi();
1499 Double_t etaTrig = ph->Eta();
1500
1501 Int_t n=fTrackEvent->GetEntriesFast() ;
1502 for(Int_t itr=0; itr<n; itr++){
0b9a6a23 1503 AliCaloPhoton * track = (AliCaloPhoton*)fTrackEvent->At(itr) ;
213dd723 1504
1505 Double_t deleta = etaTrig - track->Eta() ;
1506 Double_t delphi = phiTrig - track->Phi() ;
1507 while(delphi<-TMath::Pi()) delphi+=TMath::TwoPi() ;
1508 while(delphi>TMath::Pi()) delphi-=TMath::TwoPi() ;
1509
1510 Double_t dr = TMath::Sqrt(deleta * deleta + delphi * delphi);
1511
1512 if(dr<coneR3){
1513 eCone3+=track->Pt() ;
1514 if(dr<coneR2){
1515 eCone2+=track->Pt() ;
1516 if(dr<coneR1){
1517 eCone1+=track->Pt() ;
1518 }
1519 }
1520 }
1521 }
1522
1523 //Fill QA histgams
1524 Double_t ptTrig=ph->Pt() ;
1525 FillHistogram(Form("QA_Cone1_Tracks_cent%d",fCentBin),ptTrig,eCone1) ;
1526 FillHistogram(Form("QA_Cone2_Tracks_cent%d",fCentBin),ptTrig,eCone2) ;
1527 FillHistogram(Form("QA_Cone3_Tracks_cent%d",fCentBin),ptTrig,eCone3) ;
1528
1529
1530 //Fill Bits
1531 Int_t iCone1E1 = (epsilon1*ptTrig > eCone1) ;
1532 Int_t iCone2E1 = (epsilon1*ptTrig > eCone2) ;
1533 Int_t iCone3E1 = (epsilon1*ptTrig > eCone3) ;
1534
1535 Int_t iCone1E2 = (epsilon2*ptTrig > eCone1) ;
1536 Int_t iCone2E2 = (epsilon2*ptTrig > eCone2) ;
1537 Int_t iCone3E2 = (epsilon2*ptTrig > eCone3) ;
1538
1539
1540 Int_t isolation= iCone1E1+ 2*iCone2E1 +4*iCone3E1+
1541 8*iCone1E2 +16*iCone2E2 +32*iCone3E2 ;
1542 return isolation ;
1543}
1544//_________________________________________________________________________________
0b9a6a23 1545Bool_t AliAnalysisTaskTaggedPhotons::TestPID(Int_t iPID, AliCaloPhoton* part){
213dd723 1546// //Checks PID of particle
1547
1548 if(!part) return kFALSE ;
1549 if(iPID==0) return kTRUE ;
0b9a6a23 1550 if(iPID==1) return part->IsDispOK() ;
1551 if(iPID==2) return part->IsCPVOK() ;
1552 if(iPID==3) return part->IsDispOK() && part->IsCPVOK() ;
213dd723 1553 return kFALSE ;
1554
1555}
1556//_________________________________________________________________________________
b2a1268b 1557Double_t AliAnalysisTaskTaggedPhotons::PrimaryParticleWeight(AliAODMCParticle * particle){
1558 if(!fIsMC) //This is real data
1559 return 1;
1560 //Classify parent at vertex
1561 //Introduce for eta and pi0 weights
1562
1563 Double_t r2=particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv() ;
1564 Int_t mother = particle->GetMother() ;
1565 while(mother>-1){
1566 if(r2<1.)
1567 break ;
1568 particle = (AliAODMCParticle*) fStack->At(mother);
1569 mother = particle->GetMother() ;
1570 r2=particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv() ;
1571 }
1572
1573 //Particle within 1 cm from the virtex
1574 Int_t pdg = particle->GetPdgCode() ;
1575 Double_t x = particle->Pt() ;
1576 mother = particle->GetMother() ;
1577 while(TMath::Abs(pdg)<100){//gamma, electrons, muons
1578 if(mother>-1){
1579 AliAODMCParticle * tmpP=(AliAODMCParticle*)fStack->At(mother) ;
1580 pdg=tmpP->GetPdgCode() ;
1581 x = tmpP->Pt() ;
1582 mother = tmpP->GetMother() ;
1583 }
1584 else{ //direct photon/electron....
1585 return 1.;
1586 }
1587 }
1588 if(pdg == 111){
1589 //Pi0
1590 if(x<1) return 1. ;
1591 else return fWeightParamPi0[0]*TMath::Power(x,fWeightParamPi0[1])*
1592 (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
1593 (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) ;
1594 }
1595 return 1. ;
1596}
1597//_________________________________________________________________________________
1598void AliAnalysisTaskTaggedPhotons::SetPi0WeightParameters(TArrayD * ar){
1599
1600 for(Int_t i=0; i<6; i++){ //Array range
1601 if(ar->GetSize()>i) fWeightParamPi0[i]=ar->At(i) ;
1602 else fWeightParamPi0[i]=0.;
1603 }
1604 //normalize to obtain smooth transition at 1 GeV
1605 Double_t x=1. ;
1606 fWeightParamPi0[0]=1./(TMath::Power(x,fWeightParamPi0[1])*
1607 (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
1608 (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x)) ;
1609
213dd723 1610
1611}
1612//___________________________________________________________________________
b5f58f9c 1613Int_t AliAnalysisTaskTaggedPhotons::FindPrimary(AliVCluster*clu, Bool_t&sure){
213dd723 1614 //Finds primary and estimates if it unique one?
1615 //First check can it be photon/electron
1616 const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
1617 Int_t n=clu->GetNLabels() ;
b5f58f9c 1618 FillHistogram(Form("LabelsNPrim_cent%d",fCentBin),clu->E(),n) ;
213dd723 1619 if(n==1){
1620 sure=kTRUE;
1621 return clu->GetLabelAt(0) ;
1622 }
1623 //Number of clusters with one or more photons
1624 Bool_t hasGamma=kFALSE ;
1625 Double_t eMax=0. ;
1626 for(Int_t i=0; i<n; i++){
92071be9 1627 AliAODMCParticle* p= (AliAODMCParticle*)fStack->At(clu->GetLabelAt(i)) ;
213dd723 1628 Int_t pdg = p->GetPdgCode() ;
1629 if(pdg==22){
1630 hasGamma=kTRUE ;
92071be9 1631 if(p->E()>eMax){
1632 eMax=p->E();
213dd723 1633 }
1634 }
1635 }
1636 if(hasGamma){
b5f58f9c 1637 FillHistogram(Form("LabelsGamma_cent%d",fCentBin),clu->E()) ;
1638 FillHistogram(Form("LabelsGammaE_cent%d",fCentBin),clu->E(),eMax/clu->E()) ;
213dd723 1639 }
1640
1641 for(Int_t i=0; i<n; i++){
92071be9 1642 AliAODMCParticle* p= (AliAODMCParticle*) fStack->At(clu->GetLabelAt(i)) ;
213dd723 1643 Int_t pdg = p->GetPdgCode() ;
1644 if(pdg==22 || pdg==11 || pdg == -11){
92071be9 1645 if(p->E()>emFraction*clu->E()){
213dd723 1646 sure=kTRUE ;
1647 return clu->GetLabelAt(i);
1648 }
1649 }
1650 }
1651
1652 Double_t* Ekin= new Double_t[n] ;
1653 for(Int_t i=0; i<n; i++){
92071be9 1654 AliAODMCParticle* p=(AliAODMCParticle*) fStack->At(clu->GetLabelAt(i)) ;
213dd723 1655 Ekin[i]=p->P() ; // estimate of kinetic energy
1656 if(p->GetPdgCode()==-2212 || p->GetPdgCode()==-2112){
1657 Ekin[i]+=1.8 ; //due to annihilation
1658 }
1659 }
1660 Int_t iMax=0;
1661 Double_t eSubMax=0. ;
1662 eMax=0.;
1663 for(Int_t i=0; i<n; i++){
1664 if(Ekin[i]>eMax){
1665 eSubMax=eMax;
1666 eMax=Ekin[i];
1667 iMax=i;
1668 }
1669 }
1670 if(eSubMax>0.8*eMax)//not obvious primary
1671 sure=kFALSE;
1672 else
1673 sure=kTRUE;
1674 delete[] Ekin;
1675 return clu->GetLabelAt(iMax) ;
1676}
92071be9 1677//___________________________________________________________________________
1678Bool_t AliAnalysisTaskTaggedPhotons::IsGoodChannel(Int_t mod, Int_t ix, Int_t iz)
1679{
1680 //Check if this channel belogs to the good ones
1681
1682 if(mod>4 || mod<1){
1683 return kTRUE ;
1684 }
1685 if(!fPHOSBadMap[mod]){
1686 return kTRUE ;
1687 }
1688 if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
1689 return kFALSE ;
1690 else
1691 return kTRUE ;
1692}
213dd723 1693
1694