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