]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnalysisTaskTaggedPhotons.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
remove double initialization of emcal geometry, use virtual event and cluster classes
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnalysisTaskTaggedPhotons.cxx
CommitLineData
8b5f17b6 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:
ea7b6f24 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
d8a7052b 26// and check of the proportion of truly tagged photons.
8b5f17b6 27//
28//
29//*-- Dmitry Blau
30//////////////////////////////////////////////////////////////////////////////
31
32#include <TH1.h>
33#include <TH2.h>
34#include <TH3.h>
35#include <TFile.h>
36#include <TROOT.h>
37
38#include "AliAnalysisTaskTaggedPhotons.h"
39#include "AliAnalysisManager.h"
8b5f17b6 40#include "AliESDEvent.h"
41#include "AliAODEvent.h"
3f81d2e7 42#include "AliVCluster.h"
8b5f17b6 43#include "AliAODPWG4Particle.h"
44#include "AliAnalysisManager.h"
45#include "AliLog.h"
46#include "TGeoManager.h"
477a0971 47#include "AliMCAnalysisUtils.h"
8b5f17b6 48#include "AliMCEventHandler.h"
49#include "AliMCEvent.h"
50#include "AliStack.h"
51#include "AliPHOSGeoUtils.h"
52#include "AliEMCALGeoUtils.h"
53
b83cc5a7 54
55
8b5f17b6 56//______________________________________________________________________________
477a0971 57AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE(),
b83cc5a7 58 fPHOSgeom(0x0),
59 fEMCALgeom(0x0),
60 fStack(0x0),fPHOS(1),
477a0971 61 fPhotonId(1.0),fMinEnergyCut(0.4),
b83cc5a7 62 fPi0MeanP0(0.136),fPi0MeanP1(0.0),fPi0MeanP2(0.0),fPi0MeanP3(0.0),
ea7b6f24 63 fPi0SigmaP0(0.004508),fPi0SigmaP1(0.005497),fPi0SigmaP2(0.00000006),
477a0971 64 fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
65 fOutputList(0x0),fEventList(0x0),
38071b0f 66 fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
477a0971 67 fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
68 fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
69 fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
70 fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
b83cc5a7 71 fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
72 fhDecWMissedPartnerGeom0(0x0),
477a0971 73 fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
74 fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
75 fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
76 fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
38071b0f 77 fhMCMissedTaggingMass(0x0),
477a0971 78 fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
8b5f17b6 79{
80 //Default constructor
44b8a212 81 for (Int_t i = 0; i < 4; i++)
82 {
83 fhRecAll[i] = 0;
84 fhRecPhotonPID[i] = 0;
85 fhRecOtherPID[i] = 0;
86 fhTaggedPID[i] = 0;
87 fhInvMassReal[i] = 0;
88 fhInvMassMixed[i] = 0;
89 }
90
8b5f17b6 91}
92//______________________________________________________________________________
93AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) :
94 AliAnalysisTaskSE(name),
b83cc5a7 95 fPHOSgeom(0x0),
96 fEMCALgeom(0x0),
97 fStack(0x0),fPHOS(1),
477a0971 98 fPhotonId(1.0),fMinEnergyCut(0.4),
b83cc5a7 99 fPi0MeanP0(0.136),fPi0MeanP1(0.0),fPi0MeanP2(0.0),fPi0MeanP3(0.0),
ea7b6f24 100 fPi0SigmaP0(0.004508),fPi0SigmaP1(0.005497),fPi0SigmaP2(0.00000006),
477a0971 101 fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
102 fOutputList(0x0),fEventList(0x0),
38071b0f 103 fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
477a0971 104 fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
105 fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
106 fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
107 fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
b83cc5a7 108 fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
109 fhDecWMissedPartnerGeom0(0x0),
477a0971 110 fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
111 fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
112 fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
113 fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
38071b0f 114 fhMCMissedTaggingMass(0x0),
477a0971 115 fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
8b5f17b6 116{
117 // Constructor.
44b8a212 118 for (Int_t i = 0; i < 4; i++)
119 {
120 fhRecAll[i] = 0;
121 fhRecPhotonPID[i] = 0;
122 fhRecOtherPID[i] = 0;
123 fhTaggedPID[i] = 0;
124 fhInvMassReal[i] = 0;
125 fhInvMassMixed[i] = 0;
126 }
8b5f17b6 127
128 // Output slots
129 DefineOutput(1, TList::Class()) ;
130}
131
132//____________________________________________________________________________
133AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :
134 AliAnalysisTaskSE(ap.GetName()),
b83cc5a7 135 fPHOSgeom(0x0),
136 fEMCALgeom(0x0),
137 fStack(0x0),fPHOS(ap.fPHOS),
477a0971 138 fPhotonId(ap.fPhotonId),fMinEnergyCut(ap.fMinEnergyCut),
ea7b6f24 139 fPi0MeanP0(ap.fPi0MeanP0),fPi0MeanP1(ap.fPi0MeanP1),fPi0MeanP2(ap.fPi0MeanP2),fPi0MeanP3(ap.fPi0MeanP3),
140 fPi0SigmaP0(ap.fPi0SigmaP0),fPi0SigmaP1(ap.fPi0SigmaP1),fPi0SigmaP2(ap.fPi0SigmaP2),
477a0971 141 fZmax(ap.fZmax),fZmin(ap.fZmin),fPhimax(ap.fPhimax),fPhimin(ap.fPhimin),
142 fOutputList(0x0),fEventList(0x0),
38071b0f 143 fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
477a0971 144 fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
145 fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
146 fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
147 fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
b83cc5a7 148 fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
149 fhDecWMissedPartnerGeom0(0x0),
477a0971 150 fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
b83cc5a7 151 fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),
152 fhPartnerMissedGeo(0x0),
477a0971 153 fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
154 fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
38071b0f 155 fhMCMissedTaggingMass(0x0),
477a0971 156 fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
157{
8b5f17b6 158 // cpy ctor
44b8a212 159 for (Int_t i = 0; i < 4; i++)
160 {
161 fhRecAll[i] = 0;
162 fhRecPhotonPID[i] = 0;
163 fhRecOtherPID[i] = 0;
164 fhTaggedPID[i] = 0;
165 fhInvMassReal[i] = 0;
166 fhInvMassMixed[i] = 0;
167 }
168
8b5f17b6 169}
170
171//_____________________________________________________________________________
172AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)
173{
174// assignment operator
175
176 this->~AliAnalysisTaskTaggedPhotons();
177 new(this) AliAnalysisTaskTaggedPhotons(ap);
178 return *this;
179}
180
181//______________________________________________________________________________
182AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()
183{
184 // dtor
185 if(fOutputList) {
186 fOutputList->Clear() ;
187 delete fOutputList ;
188 }
189}
190
191
192//________________________________________________________________________
193void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
194{
195
196
197 //Load geometry
198 //if file "geometry.root" exists, load it
199 //otherwise use misaligned matrixes stored in ESD
8b5f17b6 200
201 // Create the outputs containers
202
203 OpenFile(1) ;
204 const Int_t nPtBins=52 ;
205 Double_t ptBins[nPtBins+1] ;
206 for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV: 0.1 GeV/bin
207 for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV: 0.2 GeV/bin
208 for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV: 0.5 GeV/bin
209 for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin
210 for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin
211
212
213 //Reconstructed spectra
38071b0f 214 fhRecAll[0] = new TH1D("fhRecAll0","Spectrum of all reconstructed particles, no PID", nPtBins, ptBins ) ;
215 fhRecAll[1] = new TH1D("fhRecAll1","Spectrum of all reconstructed particles, PID=1", nPtBins, ptBins ) ;
216 fhRecAll[2] = new TH1D("fhRecAll2","Spectrum of all reconstructed particles, PID=2", nPtBins, ptBins ) ;
217 fhRecAll[3] = new TH1D("fhRecAll3","Spectrum of all reconstructed particles, PID=3", nPtBins, ptBins ) ;
218
8b5f17b6 219 fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;
220 fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;
221 fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;
222
223 //Sort registered particles spectra according MC information
224 fhRecPhoton = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;
225 fhRecOther = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;
226 fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;
227 fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;
228 fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;
229 fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;
230 fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;
231 fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;
232 fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;
233 fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;
234 fhRecPhotPi0 = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;
235 fhRecPhotEta = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;
236 fhRecPhotOmega = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;
237 fhRecPhotEtapr = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;
238 fhRecPhotConv = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;
239 fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;
240 fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;
241 fhRecPhotOther = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;
242
243 //MC tagging: reasons of partner loss etc.
244 fhDecWMCPartner = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;
245 fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;
246 fhDecWMissedPartnerAll = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;
247 fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;
248 fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;
249 fhDecWMissedPartnerStack= new TH1D("fhDecWMissedPartnerStack","Decay photons with partner not in Stack", nPtBins, ptBins ) ;
250 fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;
251 fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;
252 fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;
253 fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;
254
255 //MC tagging: Decay partners spectra
256 fhPartnerMCReg = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;
257 fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;
258 fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;
259 fhPartnerMissedGeo = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;
260
261 //Tagging
262 fhTaggedAll = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;
263 fhTaggedArea1 = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;
264 fhTaggedArea2 = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;
265 fhTaggedArea3 = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;
266 fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;
267 fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;
268 fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;
269 fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;
270 fhTaggedMult = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;
271
272 //Tagging: use MC information if available
273 fhTaggedMCTrue = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;
274 fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;
275 fhMCFakeTagged = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;
276
277 //Invariant mass distributions for fake corrections
38071b0f 278 const Int_t nmass=500 ;
8b5f17b6 279 Double_t masses[nmass+1] ;
38071b0f 280 for(Int_t i=0;i<=nmass;i++)masses[i]=0.002*i ;
281 fhInvMassReal[0] = new TH2D("hInvMassRealPID0","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
282 fhInvMassReal[1] = new TH2D("hInvMassRealPID1","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
283 fhInvMassReal[2] = new TH2D("hInvMassRealPID2","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
284 fhInvMassReal[3] = new TH2D("hInvMassRealPID3","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
285 fhInvMassMixed[0] = new TH2D("hInvMassMixedPID0","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
286 fhInvMassMixed[1] = new TH2D("hInvMassMixedPID1","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
287 fhInvMassMixed[2] = new TH2D("hInvMassMixedPID2","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
288 fhInvMassMixed[3] = new TH2D("hInvMassMixedPID3","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
289 fhMCMissedTaggingMass= new TH2D("hMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;
8b5f17b6 290
291 //Conversion and annihilation radius distributions
d8a7052b 292 fhConversionRadius = new TH1D("fhConversionRadius","Radii of photon production (conversion)",100,0.,500.) ;
293 fhInteractionRadius = new TH1D("fhInteractionRadius","Radii of photon production (hadron interaction)",100,0.,500.);
8b5f17b6 294
295 fhEvents = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);
296
297// create output container
298
299 fOutputList = new TList() ;
300 fEventList = new TList() ;
301 fOutputList->SetName(GetName()) ;
302
38071b0f 303 fOutputList->AddAt(fhRecAll[0], 0) ;
304 fOutputList->AddAt(fhRecAll[1], 1) ;
305 fOutputList->AddAt(fhRecAll[2], 2) ;
306 fOutputList->AddAt(fhRecAll[3], 3) ;
307 fOutputList->AddAt(fhRecAllArea1, 4) ;
308 fOutputList->AddAt(fhRecAllArea2, 5) ;
309 fOutputList->AddAt(fhRecAllArea3, 6) ;
310
311 fOutputList->AddAt(fhRecPhoton, 7) ;
312 fOutputList->AddAt(fhRecOther, 8) ;
313 fOutputList->AddAt(fhRecPhotonPID[0], 9) ;
314 fOutputList->AddAt(fhRecPhotonPID[1], 10) ;
315 fOutputList->AddAt(fhRecPhotonPID[2], 11) ;
316 fOutputList->AddAt(fhRecPhotonPID[3], 12) ;
317 fOutputList->AddAt(fhRecOtherPID[0], 13) ;
318 fOutputList->AddAt(fhRecOtherPID[1], 14) ;
319 fOutputList->AddAt(fhRecOtherPID[2], 15) ;
320 fOutputList->AddAt(fhRecOtherPID[3], 16) ;
321 fOutputList->AddAt(fhRecPhotPi0, 17) ;
322 fOutputList->AddAt(fhRecPhotEta, 18) ;
323 fOutputList->AddAt(fhRecPhotOmega, 19) ;
324 fOutputList->AddAt(fhRecPhotEtapr, 20) ;
325 fOutputList->AddAt(fhRecPhotConv, 21) ;
326 fOutputList->AddAt(fhRecPhotHadron, 22) ;
327 fOutputList->AddAt(fhRecPhotDirect, 23) ;
328 fOutputList->AddAt(fhRecPhotOther, 24) ;
329
330 fOutputList->AddAt(fhDecWMCPartner, 25) ;
331 fOutputList->AddAt(fhDecWMissedPartnerNotPhoton, 26) ;
332 fOutputList->AddAt(fhDecWMissedPartnerAll, 27) ;
333 fOutputList->AddAt(fhDecWMissedPartnerEmin, 28) ;
334 fOutputList->AddAt(fhDecWMissedPartnerConv, 29) ;
335 fOutputList->AddAt(fhDecWMissedPartnerStack, 30) ;
336 fOutputList->AddAt(fhDecWMissedPartnerGeom0, 31) ;
337 fOutputList->AddAt(fhDecWMissedPartnerGeom1, 32) ;
338 fOutputList->AddAt(fhDecWMissedPartnerGeom2, 33) ;
339 fOutputList->AddAt(fhDecWMissedPartnerGeom3, 34) ;
340
341 fOutputList->AddAt(fhPartnerMCReg, 35) ;
342 fOutputList->AddAt(fhPartnerMissedEmin, 36) ;
343 fOutputList->AddAt(fhPartnerMissedConv, 37) ;
344 fOutputList->AddAt(fhPartnerMissedGeo, 38) ;
345
346 fOutputList->AddAt(fhTaggedAll, 39) ;
347 fOutputList->AddAt(fhTaggedArea1, 40) ;
348 fOutputList->AddAt(fhTaggedArea2, 41) ;
349 fOutputList->AddAt(fhTaggedArea3, 42) ;
350 fOutputList->AddAt(fhTaggedPID[0], 43) ;
351 fOutputList->AddAt(fhTaggedPID[1], 44) ;
352 fOutputList->AddAt(fhTaggedPID[2], 45) ;
353 fOutputList->AddAt(fhTaggedPID[3], 46) ;
354 fOutputList->AddAt(fhTaggedMult, 47) ;
355
356 fOutputList->AddAt(fhTaggedMCTrue, 48) ;
357 fOutputList->AddAt(fhMCMissedTagging, 49) ;
358 fOutputList->AddAt(fhMCFakeTagged, 50) ;
359
360 fOutputList->AddAt(fhInvMassReal[0], 51) ;
361 fOutputList->AddAt(fhInvMassReal[1], 52) ;
362 fOutputList->AddAt(fhInvMassReal[2], 53) ;
363 fOutputList->AddAt(fhInvMassReal[3], 54) ;
364 fOutputList->AddAt(fhInvMassMixed[0], 55) ;
365 fOutputList->AddAt(fhInvMassMixed[1], 56) ;
366 fOutputList->AddAt(fhInvMassMixed[2], 57) ;
367 fOutputList->AddAt(fhInvMassMixed[3], 58) ;
368 fOutputList->AddAt(fhMCMissedTaggingMass, 59) ;
369
370 fOutputList->AddAt(fhConversionRadius, 60) ;
371 fOutputList->AddAt(fhInteractionRadius, 61) ;
372
373 fOutputList->AddAt(fhEvents, 62) ;
8b5f17b6 374
375}
376
377//______________________________________________________________________________
378void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
379{
b83cc5a7 380 //Fill all histograms
c3970977 381
382
38071b0f 383 // We prepared 4 PID histograms, choose which PID combinations to use
384 // see AliAODPWG4Particle for more choises
385 // 0=all,2=Disp,4=Ch,6=Ch+Disp
386 const Int_t pidMap[4]={0,2,4,6} ;
c3970977 387
8b5f17b6 388 // Processing of one event
389 if(fDebug>1)
390 AliInfo(Form("\n\n Processing event # %lld", Entry())) ;
c3970977 391
3f81d2e7 392 AliVEvent* event = InputEvent();
393 if(!event){
394 AliDebug(1,"No event") ;
395 return;
396 }
c3970977 397
38071b0f 398 //read geometry if not read yet
399 if((fPHOS && fPHOSgeom==0) || (!fPHOS && fEMCALgeom==0))
400 InitGeometry() ;
3f81d2e7 401
402 //MC stack init
403 fStack=0x0 ;
404 if(AliAnalysisManager::GetAnalysisManager()){
405 AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
406 if(mctruth)
407 fStack = mctruth->MCEvent()->Stack();
408 }
409
410 if(!fStack && gDebug>1)
411 AliInfo("No stack! \n");
412
413
414 // Here one can choose trigger. But according to framework it should be chosen
415 // by external filters???
416 // TString trigClasses = esd->GetFiredTriggerClasses();
417 // if (!fStack && !trigClasses.Contains("CINT1B-ABCE-NOPF-ALL") ){
418 // Printf("Skip event with trigger class \"%s\"",trigClasses.Data());
419 // return;
420 // }
421
422 fhEvents->Fill(0.);
423
424 //************************ PHOS/EMCAL *************************************
425 TRefArray * caloClustersArr = new TRefArray();
426 if(fPHOS)
427 event->GetPHOSClusters(caloClustersArr);
428 else
429 event->GetEMCALClusters(caloClustersArr);
430 const Int_t kNumberOfClusters = caloClustersArr->GetEntries() ;
431
432 TClonesArray * fCaloPhotonsArr = new TClonesArray("AliAODPWG4Particle",kNumberOfClusters);
433 Int_t inList = 0; //counter of caloClusters
434
435 Int_t cluster ;
436
437 // loop over Clusters
438 for(cluster = 0 ; cluster < kNumberOfClusters ; cluster++) {
439 AliVCluster * caloCluster = static_cast<AliVCluster*>(caloClustersArr->At(cluster)) ;
c3970977 440
3f81d2e7 441 if((fPHOS && !caloCluster->IsPHOS()) ||
442 (!fPHOS && caloCluster->IsPHOS()))
443 continue ;
c3970977 444
3f81d2e7 445 Double_t v[3] ; //vertex ;
446 event->GetPrimaryVertex()->GetXYZ(v) ;
c3970977 447
3f81d2e7 448 TLorentzVector momentum ;
449 caloCluster->GetMomentum(momentum, v);
450 new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );
451 AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));
452 inList++;
c3970977 453
3f81d2e7 454 p->SetCaloLabel(cluster,-1); //This and partner cluster
455 p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));
c3970977 456
3f81d2e7 457 p->SetTag(AliMCAnalysisUtils::kMCUnknown);
458 p->SetTagged(kFALSE); //Reconstructed pairs found
459 p->SetLabel(caloCluster->GetLabel());
460 Float_t pos[3] ;
461 caloCluster->GetPosition(pos) ;
462 p->SetFiducialArea(GetFiducialArea(pos)) ;
c3970977 463
3f81d2e7 464 //PID criteria
465 p->SetDispBit(TestDisp(caloCluster->GetM02(),caloCluster->GetM20(),caloCluster->E())) ;
466 p->SetTOFBit(TestTOF(caloCluster->GetTOF(),caloCluster->E())) ;
467 p->SetChargedBit(TestCharged(caloCluster->GetEmcCpvDistance(),caloCluster->E())) ;
c3970977 468
3f81d2e7 469 for(Int_t iPID=0; iPID<4; iPID++){
470 if(p->IsPIDOK(pidMap[iPID],22))
471 fhRecAll[iPID]->Fill( p->Pt() ) ; //All recontructed particles
472 }
473 Int_t iFidArea = p->GetFiducialArea();
474 if(iFidArea>0){
475 fhRecAllArea1->Fill(p->Pt() ) ;
476 if(iFidArea>1){
477 fhRecAllArea2->Fill(p->Pt() ) ;
478 if(iFidArea>2){
479 fhRecAllArea3->Fill(p->Pt() ) ;
480 }
481 }
482 }
c3970977 483
3f81d2e7 484 if(fStack){
485 TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;
486 if(fDebug>2)
487 printf("Pdgcode = %d\n",prim->GetPdgCode());
c3970977 488
3f81d2e7 489 if(prim->GetPdgCode()!=22){ //not photon
490 //<--DP p->SetPhoton(kFALSE);
491 fhRecOther->Fill(p->Pt()); //not photons real spectra
492 for(Int_t iPID=0; iPID<4; iPID++){
493 if(p->IsPIDOK(pidMap[iPID],22))
494 fhRecOtherPID[iPID]->Fill(p->Pt());
8b5f17b6 495 }
c3970977 496 }
3f81d2e7 497 else{ //Primary photon (as in MC)
498 //<--DP p->SetPhoton(kTRUE);
499 fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon
500 for(Int_t iPID=0; iPID<4; iPID++){
501 if(p->IsPIDOK(pidMap[iPID],22))
502 fhRecPhotonPID[iPID]->Fill(p->Pt());
8b5f17b6 503 }
3f81d2e7 504 Int_t pi0i=prim->GetFirstMother();
505 Int_t grandpaPDG=-1 ;
506 TParticle * pi0p = 0;
507 if(pi0i>=0){
508 pi0p=fStack->Particle(pi0i);
509 grandpaPDG=pi0p->GetPdgCode() ;
510 }
511 switch(grandpaPDG){
512 case 111: //Pi0 decay
513 //Primary decay photon (as in MC)
514 fhRecPhotPi0->Fill(p->Pt());
515 break ;
516 case 11:
517 case -11: //electron/positron conversion
518 //<--DP p->SetConverted(1);
519 fhRecPhotConv->Fill(p->Pt()); //Reconstructed with photon from conversion primary
520 fhConversionRadius->Fill(prim->R());
521 break ;
522 case -2212:
523 case -2112: //antineutron & antiproton conversion
524 fhRecPhotHadron->Fill(p->Pt()); //Reconstructed with photon from antibaryon annihilation
525 fhInteractionRadius->Fill(prim->R());
526 break ;
527
528 case 221: //eta decay
529 fhRecPhotEta->Fill(p->Pt());
530 break ;
c3970977 531
3f81d2e7 532 case 223: //omega meson decay
533 fhRecPhotOmega->Fill(p->Pt());
534 break ;
c3970977 535
3f81d2e7 536 case 331: //eta' decay
537 fhRecPhotEtapr->Fill(p->Pt());
538 break ;
c3970977 539
3f81d2e7 540 case -1: //direct photon or no primary
541 fhRecPhotDirect->Fill(p->Pt());
542 break ;
543
544 default:
545 fhRecPhotOther->Fill(p->Pt());
546 break ;
547 }
548
549 //Now classify pi0 decay photon
550 if(grandpaPDG==111){
551 //<--DP p->Pi0Decay(kTRUE); //Mark this photon as primary decayed
552 //<--DP p->Pi0Id(pi0i); //remember id of the parent pi0
553
554 //Now check if second (partner) photon from pi0 decay hits PHOS or not
555 //i.e. both photons can be tagged or it's the systematic error
556 //<--DP p->SetPartnerPt(0.);
557 Int_t indexdecay1,indexdecay2;
558
559 indexdecay1=pi0p->GetFirstDaughter();
560 indexdecay2=pi0p->GetLastDaughter();
561 Int_t indexdecay=-1;
562 if(fDebug>2)
563 printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());
564
565 if(indexdecay1!=caloCluster->GetLabel())
566 indexdecay=indexdecay1;
567 if(indexdecay2!=caloCluster->GetLabel())
568 indexdecay=indexdecay2;
569 if(indexdecay==-1){
570 if(fDebug>2){
571 printf("Probably the other photon is not in the stack!\n");
572 printf("Number of daughters: %d\n",pi0p->GetNDaughters());
573 }
574 fhDecWMissedPartnerStack->Fill(p->Pt()) ;
575 }
576 else{
577 TParticle *partner = fStack->Particle(indexdecay);
578 //<--DP p->SetPartnerPt(partner->Pt());
579 if(partner->GetPdgCode()==22){
580 Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
581 if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector
582 if(fDebug>2)
583 printf("P_Conv, daughters=%d\n",partner->GetNDaughters());
584 //<--DP p->SetConvertedPartner(1);
585 fhPartnerMissedConv->Fill(partner->Pt());
586 fhDecWMissedPartnerConv->Fill(p->Pt()) ; //Spectrum of tagged with missed partner
587 isPartnerLost=kTRUE;
8b5f17b6 588 }
3f81d2e7 589 Bool_t impact = kFALSE ;
590 if(fPHOS){
591 Int_t modulenum ;
592 Double_t ztmp,xtmp ;
593 impact=fPHOSgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp) ;
594 if(fDebug>2){
595 printf("Impact on PHOS: module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);
b83cc5a7 596 }
3f81d2e7 597 }
598 else{
599 impact = fEMCALgeom->Impact(partner) ;
600 }
601 if(!impact){ //this photon cannot hit PHOS
602 if(fDebug>2)
603 printf("P_Geo\n");
604 fhPartnerMissedGeo->Fill(partner->Pt());
605 fhDecWMissedPartnerGeom0->Fill(p->Pt()) ; //Spectrum of tagged with missed partner
606 if(iFidArea>0){
607 fhDecWMissedPartnerGeom1->Fill(p->Pt()) ; //Spectrum of tagged with missed partner
608 if(iFidArea>1){
609 fhDecWMissedPartnerGeom2->Fill(p->Pt()) ; //Spectrum of tagged with missed partner
610 if(iFidArea>2){
611 fhDecWMissedPartnerGeom3->Fill(p->Pt()) ; //Spectrum of tagged with missed partner
8b5f17b6 612 }
613 }
614 }
3f81d2e7 615 isPartnerLost=kTRUE;
616 }
617 if(!isPartnerLost && partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS
618 if(fDebug>2)
619 printf("P_Reg, E=%f\n",partner->Energy());
620 fhPartnerMissedEmin->Fill(partner->Pt()); //Spectrum of missed partners
621 fhDecWMissedPartnerEmin->Fill(p->Pt()) ; //Spectrum of tagged with missed partner
622 isPartnerLost=kTRUE;
623 }
624 if(!isPartnerLost){
625 // p->SetMCTagged(1); //set this photon as primary tagged
626 fhDecWMCPartner->Fill(p->Pt());
627 fhPartnerMCReg->Fill(partner->Pt());
628 if(fDebug>2){
629 printf("both photons are inside PHOS. Energy: %f, Pt of pair photon: %f, E of pair photon: %f, Px: %f Py: %f Pz: %f, num of daughters: %d \n", caloCluster->E(),partner->Pt(),partner->Energy(),partner->Px(),partner->Py(),partner->Pz(),partner->GetNDaughters());
c3970977 630 }
8b5f17b6 631 }
3f81d2e7 632 else{
633 fhDecWMissedPartnerAll->Fill(p->Pt());
634 }
635 }//Partner - photon
636 else{//partner not photon
637 fhDecWMissedPartnerNotPhoton->Fill(p->Pt());
8b5f17b6 638 }
639 }
640 }
641 }
3f81d2e7 642 }
643 } //PHOS/EMCAL clusters
644
645 if(fDebug>1)
646 printf("number of clusters: %d\n",inList);
647
648 //Invariant Mass analysis
649 for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList ; phosPhoton1++) {
650 AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1));
651 for(Int_t phosPhoton2 = 0 ; phosPhoton2 < inList ; phosPhoton2++) {
652 if(phosPhoton1==phosPhoton2)
653 continue ;
654 AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));
655
656 Double_t invMass = p1->GetPairMass(p2);
657 for(Int_t iPID=0; iPID<4; iPID++){
658 if(p1->IsPIDOK(pidMap[iPID],22))
659 fhInvMassReal[iPID]->Fill(invMass,p1->Pt());
660 if(p2->IsPIDOK(pidMap[iPID],22))
661 fhInvMassReal[iPID]->Fill(invMass,p2->Pt());
662 }
663 if(fDebug>2)
664 printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);
665
666 Bool_t makePi0=IsInPi0Band(invMass,p1->Pt());
667
668 if(makePi0 && p1->IsTagged()){//Multiple tagging
669 fhTaggedMult->Fill(p1->Pt());
670 }
671 if(makePi0 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times
672 fhTaggedAll->Fill(p1->Pt());
673 Int_t iFidArea = p1->GetFiducialArea();
674 if(iFidArea>0){
675 fhTaggedArea1->Fill(p1->Pt() ) ;
676 if(iFidArea>1){
677 fhTaggedArea2->Fill(p1->Pt() ) ;
678 if(iFidArea>2){
679 fhTaggedArea3->Fill(p1->Pt() ) ;
680 }
681 }
682 }
c3970977 683
c3970977 684 for(Int_t iPID=0; iPID<4; iPID++){
685 if(p1->IsPIDOK(pidMap[iPID],22))
3f81d2e7 686 fhTaggedPID[iPID]->Fill(p1->Pt());
c3970977 687 }
c3970977 688
3f81d2e7 689 p1->SetTagged(kTRUE) ;
690 }
691
692 //First chesk if this is true pi0 pair
693 if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0
694 // p1->SetTrueTagged(1);
695 // p2->SetTrueTagged(1);
696 if(makePi0)//Correctly tagged photons
697 fhTaggedMCTrue->Fill(p1->Pt());
698 else{ //Decay pair missed tagging
699 fhMCMissedTagging->Fill(p1->Pt());
700 fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;
701 //Clussify why missed tagging (todo)
702 //Converted
703 //Partner not a photon
704 //Tagged not a photon
705 //Just wrong inv.mass
c3970977 706 }
3f81d2e7 707 }
708 else{//Fake tagged - not from the same pi0
709 if(makePi0)//Fake pair
710 fhMCFakeTagged->Fill(p1->Pt());
8b5f17b6 711 }
712 }
3f81d2e7 713 }
714
715 //Fill Mixed InvMass distributions:
716 TIter nextEv(fEventList) ;
717 while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
718 Int_t nPhotons2 = event2->GetEntriesFast() ;
719 for(Int_t i=0; i < inList ; i++){
720 AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;
721 for(Int_t j=0; j < nPhotons2 ; j++){
722 AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;
723 Double_t invMass = p1->GetPairMass(p2);
724 for(Int_t iPID=0; iPID<4; iPID++){
725 if(p1->IsPIDOK(pidMap[iPID],22))
726 fhInvMassMixed[iPID]->Fill(invMass,p1->Pt());
727 if(p2->IsPIDOK(pidMap[iPID],22))
728 fhInvMassMixed[iPID]->Fill(invMass,p2->Pt());
38071b0f 729 }
8b5f17b6 730 }
731 }
3f81d2e7 732 }
733
734 //Remove old events
735 fEventList->AddFirst(fCaloPhotonsArr);
736 if(fEventList->GetSize() > 10){
737 TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());
738 fEventList->Remove(tmp);
739 delete tmp;
740 }
741
742 delete caloClustersArr;
743
744 PostData(1, fOutputList);
745
8b5f17b6 746}
747
748
3f81d2e7 749
8b5f17b6 750//______________________________________________________________________________
751void AliAnalysisTaskTaggedPhotons::Init()
752{
753 // Intialisation of parameters
754 AliInfo("Doing initialisation") ;
755 SetPhotonId(0.9) ;
756 SetMinEnergyCut(0.4);
757 SetPi0MeanParameters(0.136,0.,0.0,0.0);
758// SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);
759 SetPi0SigmaParameters(0.004508,0.005497,0.00000006);
760}
761
762//______________________________________________________________________________
763void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
764{
765 // Processing when the event loop is ended
44b8a212 766 if (fDebug > 1) Printf("Terminate()");
8b5f17b6 767}
768//______________________________________________________________________________
769Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
770{
771 //Parameterization of the pi0 peak region
772 Double_t mpi0mean = fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;
773 Double_t mpi0sigma = TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);
774
775 return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
776}
777//______________________________________________________________________________
778Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{
779 //Looks through parents and finds if there was commont pi0 among ancestors
780
781 if(!fStack)
782 return kFALSE ; //can not say anything
783
784 Int_t prim1 = p1->GetLabel();
785 while(prim1!=-1){
786 Int_t prim2 = p2->GetLabel();
787 while(prim2!=-1){
788 if(prim1==prim2){
789 if(fStack->Particle(prim1)->GetPdgCode()==111)
790 return kTRUE ;
791 else
792 return kFALSE ;
793 }
794 prim2=fStack->Particle(prim2)->GetFirstMother() ;
795 }
796 prim1=fStack->Particle(prim1)->GetFirstMother() ;
797 }
798 return kFALSE ;
799}
800//______________________________________________________________________________
801Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos)const{
802 //calculates in which kind of fiducial area photon hit
803 Double_t phi=TMath::ATan2(pos[1],pos[0]) ;
b83cc5a7 804 Double_t z=pos[2] ;
8b5f17b6 805 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
806 while(phi<0.)phi+=TMath::TwoPi() ;
917d08eb 807 if(fDebug>2){
808 printf("FiducialArea: phi=%f, z=%f \n",phi,z) ;
809 printf(" fZmax=%f, fZmin=%f, fPhimax=%f, fPhimin=%f \n",fZmax,fZmin,fPhimax,fPhimin) ;
810 }
8b5f17b6 811 if(fPHOS){
812 //From active PHOS area remove bands in 10 cm
813 const Double_t kphi=TMath::ATan(10./460.) ; //angular band width
814 Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;
815 Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;
816 Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
817 Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
917d08eb 818 if(fDebug>2){
819 printf("In PHOS \n") ;
820 printf(" dzMax=%f, dzMin=%f, dphiMax=%f, dphiMin=%f ret=%d\n",dzMax,dzMin,dphiMax,dphiMin,(Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin))) ;
821 }
8b5f17b6 822 return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin));
823 }
824 else{//EMCAL
825 //From active EMCAL area remove bands in 20 cm
826 const Double_t kphi=TMath::ATan(20./428.) ; //angular band width
827 Double_t dzMax=TMath::Ceil((fZmax-z)/20.) ;
828 Double_t dzMin=TMath::Ceil((z-fZmin)/20.) ;
829 Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
830 Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
831 return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin));
832 }
833}
834//______________________________________________________________________________
38071b0f 835Bool_t AliAnalysisTaskTaggedPhotons::TestDisp(Double_t l0, Double_t l1, Double_t /*e*/)const{
b83cc5a7 836 //test if dispersion corresponds to those of photon
837 if(fPHOS){
3f81d2e7 838 // Double_t l0mean=1.38736+0.490405*TMath::Exp(-e*0.286170) ;
839 // Double_t l1mean=1.09786-0.323469*TMath::Exp(-e*0.918719) ;
840 // Double_t l0sigma=0.159905+0.829831/e-0.158067/e/e ;
841 // Double_t l1sigma=0.133170+0.404387/e-0.0426302/e/e ;
842 // Double_t c =-0.382233 ;
843 // return ((l0-l0mean)*(l0-l0mean)/l0sigma/l0sigma + (l1-l1mean)*(l1-l1mean)/l1sigma/l1sigma+c*(l0-l0mean)*(l1-l1mean)/l0sigma/l1sigma)<1. ;
844
845 if(l1>= 0 && l0>= 0 && l1 < 0.1 && l0 < 0.1) return kFALSE;
846 if(l1>= 0 && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) return kTRUE;
847 if(l1>= 0 && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) return kFALSE;
848 if(l1>= 0 && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) return kFALSE;
849 if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) return kTRUE;
850 if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) return kFALSE;
851 if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) return kFALSE;
852 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) return kTRUE;
853 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) return kTRUE;
854 if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) return kTRUE;
855 return kFALSE ;
856
b83cc5a7 857 }
858 else{ //EMCAL: not ready yet
3f81d2e7 859 return kTRUE ;
860
b83cc5a7 861 }
3f81d2e7 862
b83cc5a7 863}
38071b0f 864//______________________________________________________________________________^M
865Bool_t AliAnalysisTaskTaggedPhotons::TestCharged(Double_t dr,Double_t /*en*/)const{
3f81d2e7 866 // Test charged
867
38071b0f 868 if(dr<15.) return kFALSE ;
869 return kTRUE ;
b83cc5a7 870
38071b0f 871}
872//______________________________________________________________________________^M
873void AliAnalysisTaskTaggedPhotons::InitGeometry(){
874 //Read rotation matrixes from ESD
c3970977 875
917d08eb 876 if(fDebug>1)printf("Init geometry \n") ;
3f81d2e7 877 AliESDEvent* esd = dynamic_cast<AliESDEvent*> (InputEvent()) ;
878 AliAODEvent* aod = 0x0;
879 if(!esd) aod = dynamic_cast<AliAODEvent*> (InputEvent()) ;
880
881 if(!aod && !esd) return;
882
883 if(fPHOS){//reading PHOS matrixes
884 fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
885 Bool_t activeMod[5]={0} ;
886 for(Int_t mod=0; mod<5; mod++){
887 if(esd){
888 const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
889 fPHOSgeom->SetMisalMatrix(m, mod) ;
890 if(m!=0)
891 activeMod[mod]=kTRUE ;
38071b0f 892 }
3f81d2e7 893 else{
894 const TGeoHMatrix* m=aod->GetHeader()->GetPHOSMatrix(mod) ;
895 fPHOSgeom->SetMisalMatrix(m, mod) ;
896 if(m!=0)
897 activeMod[mod]=kTRUE ;
38071b0f 898 }
899 }
3f81d2e7 900
901 if(fZmax>fZmin) //already set manually
902 return ;
903
904 fZmax= 999. ;
905 fZmin=-999. ;
906 fPhimax=-999. ;
907 fPhimin= 999. ;
908 for(Int_t imod=1; imod<=5; imod++){
909 if(!activeMod[imod-1])
910 continue ;
911 //Find exact coordinates of PHOS corners
912 Int_t relId[4]={imod,0,1,1} ;
913 Float_t x,z ;
914 fPHOSgeom->RelPosInModule(relId,x,z) ;
915 TVector3 pos ;
916 fPHOSgeom->Local2Global(imod,x,z,pos) ;
917 Double_t phi=pos.Phi() ;
918 while(phi<0.)phi+=TMath::TwoPi() ;
919 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
920 fPhimin=TMath::Min(fPhimin,float(phi)) ;
921 fZmin=TMath::Max(fZmin,float(pos.Z())) ;
922 relId[2]=64 ;
923 relId[3]=56 ;
924 fPHOSgeom->RelPosInModule(relId,x,z) ;
925 fPHOSgeom->Local2Global(imod,x,z,pos) ;
926 phi=pos.Phi() ;
927 while(phi<0.)phi+=TMath::TwoPi() ;
928 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
929 fPhimax=TMath::Max(fPhimax,float(phi)) ;
930 fZmax=TMath::Min(fZmax,float(pos.Z())) ;
931 }
932 }
933 else{ //EMCAL
934 fEMCALgeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1");
935 for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
936 if(esd){
937 const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
938 fEMCALgeom->SetMisalMatrix(m, mod) ;
38071b0f 939 }
3f81d2e7 940 else{
941 const TGeoHMatrix* m=aod->GetHeader()->GetEMCALMatrix(mod) ;
942 fEMCALgeom->SetMisalMatrix(m, mod) ;
38071b0f 943 }
944 }
3f81d2e7 945 if(fZmax>fZmin) //already set manually
946 return ;
947
948 fZmax= 999. ;
949 fZmin=-999. ;
950 fPhimax=-999. ;
951 fPhimin= 999. ;
952 for(Int_t imod=0; imod<12; imod++){
953 //Find exact coordinates of SM corners
954 Int_t absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 0, 0);
955 TVector3 pos ;
956 //Get the position of this tower.
957 fEMCALgeom->RelPosCellInSModule(absId,pos);
958 Double_t phi=pos.Phi() ;
959 while(phi<0.)phi+=TMath::TwoPi() ;
960 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
961 fPhimin=TMath::Min(fPhimin,float(phi)) ;
962 fZmin=TMath::Max(fZmin,float(pos.Z())) ;
963 absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 24, 48);
964 fEMCALgeom->RelPosCellInSModule(absId,pos);
965 phi=pos.Phi() ;
966 while(phi<0.)phi+=TMath::TwoPi() ;
967 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
968 fPhimax=TMath::Max(fPhimax,float(phi)) ;
969 fZmax=TMath::Min(fZmax,float(pos.Z())) ;
970 }
971 }
38071b0f 972}