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