]> 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
new task for direct photon identification tagging photons as decay or not
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnalysisTaskTaggedPhotons.cxx
CommitLineData
477a0971 1/**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15\r
16/* $Id: */\r
17\r
18//_________________________________________________________________________\r
19// Analysis for Tagged Photons\r
20// Prepares all necessary histograms for calculation of \r
21// the yield of pi0 decay photon in calorimeter:\r
22// \r
23//\r
24//*-- Dmitry Blau \r
25//////////////////////////////////////////////////////////////////////////////\r
26\r
27#include <TCanvas.h>\r
28#include <TH1.h>\r
29#include <TH2.h>\r
30#include <TH3.h>\r
31#include <TROOT.h>\r
32#include <TLegend.h> \r
33#include <TNtuple.h>\r
34#include <TVector3.h> \r
35#include <TFile.h>\r
36\r
37#include "AliAnalysisTaskTaggedPhotons.h" \r
38#include "AliAnalysisManager.h"\r
39#include "AliESDEvent.h" \r
40#include "AliESDCaloCluster.h" \r
41#include "AliAODPWG4Particle.h"\r
42#include "AliLog.h"\r
43#include "AliESDVertex.h"\r
44#include "AliPHOSGeoUtils.h"\r
45#include "TGeoManager.h"\r
46#include "AliMCAnalysisUtils.h"
47#include "AliMCEventHandler.h"\r
48#include "AliAnalysisManager.h"\r
49#include "AliMCEvent.h"\r
50#include "AliStack.h"\r
51\r
52//______________________________________________________________________________\r
53AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE(),
54 fgeom(0x0),fStack(0x0),fDebug(0),fPHOS(1),
55 fPhotonId(1.0),fMinEnergyCut(0.4),
56 fPi0Mean_p0(0.1377),fPi0Mean_p1(-0.002566),fPi0Mean_p2(0.001216),fPi0Mean_p3(-0.0001256),
57 fPi0Sigma_p0(0.004508),fPi0Sigma_p1(0.005497),fPi0Sigma_p2(0.00000006),
58 fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
59 fOutputList(0x0),fEventList(0x0),
60 fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
61 fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
62 fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
63 fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
64 fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
65 fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerGeom0(0x0),
66 fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
67 fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
68 fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
69 fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
70 fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
71 fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
72{\r
73 //Default constructor\r
74}\r
75//______________________________________________________________________________\r
76AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) : \r
77 AliAnalysisTaskSE(name),\r
78 fgeom(0x0),fStack(0x0),fDebug(0),fPHOS(1),
79 fPhotonId(1.0),fMinEnergyCut(0.4),
80 fPi0Mean_p0(0.1377),fPi0Mean_p1(-0.002566),fPi0Mean_p2(0.001216),fPi0Mean_p3(-0.0001256),
81 fPi0Sigma_p0(0.004508),fPi0Sigma_p1(0.005497),fPi0Sigma_p2(0.00000006),
82 fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
83 fOutputList(0x0),fEventList(0x0),
84 fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
85 fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
86 fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
87 fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
88 fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
89 fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerGeom0(0x0),
90 fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
91 fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
92 fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
93 fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
94 fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
95 fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
96{\r
97 // Constructor.\r
98\r
99 // Output slots \r
100 DefineOutput(1, TList::Class()) ; \r
101}\r
102\r
103//____________________________________________________________________________\r
104AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :\r
105 AliAnalysisTaskSE(ap.GetName()), \r
106 fgeom(0x0),fStack(0x0),fDebug(ap.fDebug),fPHOS(ap.fPHOS),
107 fPhotonId(ap.fPhotonId),fMinEnergyCut(ap.fMinEnergyCut),
108 fPi0Mean_p0(ap.fPi0Mean_p0),fPi0Mean_p1(ap.fPi0Mean_p1),fPi0Mean_p2(ap.fPi0Mean_p2),fPi0Mean_p3(ap.fPi0Mean_p3),
109 fPi0Sigma_p0(ap.fPi0Sigma_p0),fPi0Sigma_p1(ap.fPi0Sigma_p1),fPi0Sigma_p2(ap.fPi0Sigma_p2),
110 fZmax(ap.fZmax),fZmin(ap.fZmin),fPhimax(ap.fPhimax),fPhimin(ap.fPhimin),
111 fOutputList(0x0),fEventList(0x0),
112 fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
113 fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
114 fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
115 fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
116 fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
117 fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerGeom0(0x0),
118 fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
119 fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
120 fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
121 fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
122 fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
123 fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
124{
125 // cpy ctor\r
126}\r
127\r
128//_____________________________________________________________________________\r
129AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)\r
130{\r
131// assignment operator\r
132\r
133 this->~AliAnalysisTaskTaggedPhotons();\r
134 new(this) AliAnalysisTaskTaggedPhotons(ap);\r
135 return *this;\r
136}\r
137\r
138//______________________________________________________________________________\r
139AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()\r
140{\r
141 // dtor\r
142 if(fOutputList) {\r
143 fOutputList->Clear() ; \r
144 delete fOutputList ;\r
145 }\r
146}\r
147\r
148\r
149//________________________________________________________________________\r
150void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()\r
151{ \r
152\r
153 //Load geometry\r
154 //if file "geometry.root" exists, load it\r
155 //otherwise use misaligned matrixes stored in ESD\r
156 TFile *geoFile = new TFile("geometry.root","read");\r
157 if(geoFile->IsZombie()){ //no file, load geo matrixes from ESD\r
158 //todo\r
159 AliInfo("Can not find file geometry.root, reading misalignment matrixes from AliESDs") ;
160 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
161 if(!esd)
162 AliFatal("Can not read geometry even from ESD. Note, that AOD does not contain PHOS/EMCAL geometry") ;
163 if(fPHOS){//reading PHOS matrixes
164 fgeom = new AliPHOSGeoUtils("IHEP","");\r
165 for(Int_t mod=0; mod<5; mod++){
166 const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
167 fgeom->SetMisalMatrix(m, mod) ;
168 }
169 }
170 }\r
171 else{\r
172 gGeoManager = (TGeoManager*)geoFile->Get("Geometry");\r
173 fgeom = new AliPHOSGeoUtils("IHEP","");\r
174 }\r
175\r
176\r
177 if(fgeom==NULL){\r
178 AliError("Error loading Geometry\n");\r
179 }\r
180 else\r
181 AliInfo("Geometry loaded... OK\n");\r
182\r
183 //Evaluate active PHOS/EMCAL area\r
184 //To be fixed todo\r
185 fZmax= 2.25*56/2. ;\r
186 fZmin=-2.25*56/2. ;\r
187 fPhimax=220./180.*TMath::Pi() ;\r
188 fPhimin=320./180.*TMath::Pi() ;\r
189\r
190\r
191 // Create the outputs containers\r
192\r
193 OpenFile(1) ; \r
194 const Int_t nPtBins=51 ;\r
195 Double_t ptBins[nPtBins+1] ;\r
196 for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV: 0.1 GeV/bin\r
197 for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV: 0.2 GeV/bin \r
198 for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV: 0.5 GeV/bin \r
199 for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin \r
200 for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin \r
201\r
202\r
203 //Reconstructed spectra\r
204 fhRecAll = new TH1D("fhRecAll","Spectrum of all reconstructed particles", nPtBins, ptBins ) ;\r
205 fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;\r
206 fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;\r
207 fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;\r
208\r
209 //Sort registered particles spectra according MC information\r
210 fhRecPhoton = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;\r
211 fhRecOther = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;\r
212 fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;\r
213 fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;\r
214 fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;\r
215 fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;\r
216 fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;\r
217 fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;\r
218 fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;\r
219 fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;\r
220 fhRecPhotPi0 = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;\r
221 fhRecPhotEta = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;\r
222 fhRecPhotOmega = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;\r
223 fhRecPhotEtapr = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;\r
224 fhRecPhotConv = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;\r
225 fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;\r
226 fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;\r
227 fhRecPhotOther = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;\r
228\r
229 //MC tagging: reasons of partner loss etc.\r
230 fhDecWMCPartner = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;\r
231 fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;\r
232 fhDecWMissedPartnerAll = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;\r
233 fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;\r
234 fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;\r
235 fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;\r
236 fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;\r
237 fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;\r
238 fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;\r
239\r
240 //MC tagging: Decay partners spectra\r
241 fhPartnerMCReg = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;\r
242 fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;\r
243 fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;\r
244 fhPartnerMissedGeo = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;\r
245\r
246 //Tagging\r
247 fhTaggedAll = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;\r
248 fhTaggedArea1 = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;\r
249 fhTaggedArea2 = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;\r
250 fhTaggedArea3 = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;\r
251 fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;\r
252 fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;\r
253 fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;\r
254 fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;\r
255 fhTaggedMult = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;\r
256\r
257 //Tagging: use MC information if available\r
258 fhTaggedMCTrue = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;\r
259 fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;\r
260 fhMCFakeTagged = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;\r
261\r
262 //Invariant mass distributions for fake corrections\r
263 const Int_t nmass=200 ;\r
264 Double_t masses[nmass+1] ;\r
265 for(Int_t i=0;i<=nmass;i++)masses[i]=0.005*i ;\r
266 fhInvMassReal = new TH2D("fhInvMassReal","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
267 fhInvMassMixed = new TH2D("fhInvMassMixed","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
268 fhMCMissedTaggingMass= new TH2D("fhMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;\r
269\r
270 //Conversion and annihilation radius distributions\r
271 fhConversionRadius = new TH1D("fhConversionRadius","Radis of photon production (conversion)",100,0.,500.) ;\r
272 fhInteractionRadius = new TH1D("fhInteractionRadius","Radis of photon production (hadron interaction)",100,0.,500.);\r
273\r
274 fhEvents = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);\r
275\r
276// create output container\r
277 \r
278 fOutputList = new TList() ;\r
279 fEventList = new TList() ; \r
280 fOutputList->SetName(GetName()) ; \r
281\r
282 fOutputList->AddAt(fhRecAll, 0) ;\r
283 fOutputList->AddAt(fhRecAllArea1, 1) ;\r
284 fOutputList->AddAt(fhRecAllArea2, 2) ;\r
285 fOutputList->AddAt(fhRecAllArea3, 3) ;\r
286\r
287 fOutputList->AddAt(fhRecPhoton, 4) ;\r
288 fOutputList->AddAt(fhRecOther, 5) ;\r
289 fOutputList->AddAt(fhRecPhotonPID[0], 6) ;\r
290 fOutputList->AddAt(fhRecPhotonPID[1], 7) ;\r
291 fOutputList->AddAt(fhRecPhotonPID[2], 8) ;\r
292 fOutputList->AddAt(fhRecPhotonPID[3], 9) ;\r
293 fOutputList->AddAt(fhRecOtherPID[0], 10) ;\r
294 fOutputList->AddAt(fhRecOtherPID[1], 11) ;\r
295 fOutputList->AddAt(fhRecOtherPID[2], 12) ;\r
296 fOutputList->AddAt(fhRecOtherPID[3], 13) ;\r
297 fOutputList->AddAt(fhRecPhotPi0, 14) ;\r
298 fOutputList->AddAt(fhRecPhotEta, 15) ;\r
299 fOutputList->AddAt(fhRecPhotOmega, 16) ;\r
300 fOutputList->AddAt(fhRecPhotEtapr, 17) ;\r
301 fOutputList->AddAt(fhRecPhotConv, 18) ;\r
302 fOutputList->AddAt(fhRecPhotHadron, 19) ;\r
303 fOutputList->AddAt(fhRecPhotDirect, 20) ;\r
304 fOutputList->AddAt(fhRecPhotOther, 21) ;\r
305\r
306 fOutputList->AddAt(fhDecWMCPartner, 22) ;\r
307 fOutputList->AddAt(fhDecWMissedPartnerNotPhoton, 23) ;\r
308 fOutputList->AddAt(fhDecWMissedPartnerAll, 24) ;\r
309 fOutputList->AddAt(fhDecWMissedPartnerEmin, 25) ;\r
310 fOutputList->AddAt(fhDecWMissedPartnerConv, 26) ;\r
311 fOutputList->AddAt(fhDecWMissedPartnerGeom0, 27) ;\r
312 fOutputList->AddAt(fhDecWMissedPartnerGeom1, 28) ;\r
313 fOutputList->AddAt(fhDecWMissedPartnerGeom2, 29) ;\r
314 fOutputList->AddAt(fhDecWMissedPartnerGeom3, 30) ;\r
315\r
316 fOutputList->AddAt(fhPartnerMCReg, 31) ;\r
317 fOutputList->AddAt(fhPartnerMissedEmin, 32) ;\r
318 fOutputList->AddAt(fhPartnerMissedConv, 33) ;\r
319 fOutputList->AddAt(fhPartnerMissedGeo, 34) ;\r
320\r
321 fOutputList->AddAt(fhTaggedAll, 35) ;\r
322 fOutputList->AddAt(fhTaggedArea1, 36) ;\r
323 fOutputList->AddAt(fhTaggedArea2, 37) ;\r
324 fOutputList->AddAt(fhTaggedArea3, 38) ;\r
325 fOutputList->AddAt(fhTaggedPID[0], 39) ;\r
326 fOutputList->AddAt(fhTaggedPID[1], 40) ;\r
327 fOutputList->AddAt(fhTaggedPID[2], 41) ;\r
328 fOutputList->AddAt(fhTaggedPID[3], 42) ;\r
329 fOutputList->AddAt(fhTaggedMult, 43) ;\r
330\r
331 fOutputList->AddAt(fhTaggedMCTrue, 44) ;\r
332 fOutputList->AddAt(fhMCMissedTagging, 45) ;\r
333 fOutputList->AddAt(fhMCFakeTagged, 46) ;\r
334\r
335 fOutputList->AddAt(fhInvMassReal, 47) ;\r
336 fOutputList->AddAt(fhInvMassMixed, 48) ;\r
337 fOutputList->AddAt(fhMCMissedTaggingMass, 49) ;\r
338\r
339 fOutputList->AddAt(fhConversionRadius, 50) ;\r
340 fOutputList->AddAt(fhInteractionRadius, 51) ;\r
341\r
342 fOutputList->AddAt(fhEvents, 52) ;\r
343\r
344\r
345/*\r
346 fOutputList->AddAt(fhPHOSPos, 0) ; \r
347 fOutputList->AddAt(fhPHOS, 1) ; \r
348 fOutputList->AddAt(fhAllPhotons, 2) ;\r
349 fOutputList->AddAt(fhNotPhotons, 3) ;\r
350 fOutputList->AddAt(fhAllPhotonsPrimary, 4) ;\r
351 fOutputList->AddAt(fhNotPhotonsPrimary, 5) ;\r
352 fOutputList->AddAt(fhfakeNotPhotons, 6) ;\r
353\r
354 fOutputList->AddAt(fhTaggedPhotons, 7) ; \r
355 fOutputList->AddAt(fhfakeTaggedPhotons, 8) ;\r
356 fOutputList->AddAt(fhDecayNotTaggedPhotons, 9) ;\r
357 fOutputList->AddAt(fhPi0DecayPhotonsPrimary, 10);\r
358 fOutputList->AddAt(fhEtaDecayPhotonsPrimary, 11);\r
359 fOutputList->AddAt(fhOmegaDecayPhotonsPrimary, 12);\r
360 fOutputList->AddAt(fhEtaSDecayPhotonsPrimary, 13);\r
361 fOutputList->AddAt(fhOtherDecayPhotonsPrimary, 14);\r
362 fOutputList->AddAt(fhDecayPhotonsPrimary, 15);\r
363\r
364 fOutputList->AddAt(fhConvertedPhotonsPrimary, 16);\r
365 fOutputList->AddAt(fhConvertedPhotonsPrimaryHadronsDecays, 17);\r
366 fOutputList->AddAt(fhCoordsConversion, 18);\r
367 fOutputList->AddAt(fhCoordsConversion2, 19);\r
368\r
369 fOutputList->AddAt(fhPHOSInvariantMassReal, 20);\r
370 fOutputList->AddAt(fhPHOSInvariantMassMixed, 21);\r
371\r
372 fOutputList->AddAt(fhPHOSPi0, 22);\r
373\r
374 fOutputList->AddAt(fhPi0DecayPhotonsGeomfake, 23);\r
375 fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimary, 24);\r
376 fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimaryPair, 25);\r
377 fOutputList->AddAt(fhPi0DecayPhotonsBigDecay, 26);\r
378 fOutputList->AddAt(fhPi0DecayPhotonsPConv, 27);\r
379 fOutputList->AddAt(fhPi0DecayPhotonsPGeo, 28);\r
380 fOutputList->AddAt(fhPi0DecayPhotonsPReg, 29);\r
381\r
382 fOutputList->AddAt(fhfakeTaggedPhotonsConv, 30) ;\r
383 fOutputList->AddAt(fhfakeTaggedPhotonsPID, 31) ;\r
384 fOutputList->AddAt(fhstrangeNotTaggedPhotons, 32) ;\r
385 fOutputList->AddAt(fhstrangeNotTaggedPhotonsPair, 33) ;\r
386 fOutputList->AddAt(fhstrangeNotTaggedPhotonsRegCut, 34) ;\r
387 fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimaryRegCut, 35);\r
388 fOutputList->AddAt(fhstrangeNotTaggedPhotonsPairRegCut, 36) ;\r
389 fOutputList->AddAt(fhPi0DecayPhotonsTaggedPrimaryPairRegCut, 37);\r
390 fOutputList->AddAt(fhTrackRefCoords, 38);\r
391*/\r
392 fOutputList->AddAt(fhEvents, 39);\r
393}\r
394\r
395//______________________________________________________________________________\r
396void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *) \r
397{\r
398 fhEvents->Fill(0.);\r
399\r
400 // Processing of one event\r
401 if(fDebug>1)\r
402 AliInfo(Form("\n\n Processing event # %lld", Entry())) ; \r
403 AliESDEvent* esd = (AliESDEvent*)InputEvent();\r
404\r
405 if(fDebug>2){ //DP: Check these histograms <------- \r
406// printf("Tagged: %f ",fhTaggedPhotons->GetEntries());\r
407// printf("fakeTagged: %f ",fhfakeTaggedPhotons->GetEntries());\r
408// printf("fakeNotTagged: %f ",fhfakeNotTaggedPhotons->GetEntries());\r
409// printf("strangeNotTagged: %f ",fhstrangeNotTaggedPhotons->GetEntries());\r
410// printf("trueTagged: %f ",fhPi0DecayPhotonsTrueTagged->GetEntries());\r
411// printf("fakeTaggedConv: %f ",fhfakeTaggedPhotonsConv->GetEntries());\r
412// printf("\nDIFF: %f\n",fhTaggedPhotons->GetEntries()-fhfakeTaggedPhotons->GetEntries()+fhfakeNotTaggedPhotons->GetEntries()+fhstrangeNotTaggedPhotons->GetEntries()-fhPi0DecayPhotonsTrueTagged->GetEntries()-fhfakeTaggedPhotonsConv->GetEntries());\r
413 }\r
414\r
415 //MC stack init\r
416 AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
417 fStack = mctruth->MCEvent()->Stack();\r
418\r
419 if(!fStack && gDebug>1)\r
420 AliInfo("No stack! \n");\r
421\r
422 //************************ PHOS *************************************\r
423 TRefArray * caloClustersArr = new TRefArray(); \r
424 esd->GetPHOSClusters(caloClustersArr);\r
425 const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ; \r
426\r
427 TClonesArray * fCaloPhotonsArr = new TClonesArray("AliAODPWG4Particle",kNumberOfPhosClusters);\r
428 Int_t inList = 0; //counter of caloClusters\r
429\r
430 Int_t phosCluster ; \r
431\r
432 // loop over PHOS Clusters\r
433 for(phosCluster = 0 ; phosCluster < kNumberOfPhosClusters ; phosCluster++) {\r
434 AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(phosCluster)) ;\r
435 \r
436 if((fPHOS && !caloCluster->IsPHOS()) ||\r
437 (!fPHOS && caloCluster->IsPHOS()))\r
438 continue ; \r
439\r
440 Double_t v[3] ; //vertex ;\r
441 esd->GetVertex()->GetXYZ(v) ;\r
442\r
443 TLorentzVector momentum ;\r
444 caloCluster->GetMomentum(momentum, v);\r
445 new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );\r
446 AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));\r
447 inList++;\r
448\r
449 p->SetCaloLabel(phosCluster,-1); //This and partner cluster
450 p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));\r
451
452 p->SetTag(AliMCAnalysisUtils::kMCUnknown);\r
453 p->SetTagged(kFALSE); //Reconstructed pairs found\r
454 p->SetLabel(caloCluster->GetLabel());\r
455 Float_t pos[3] ;
456 caloCluster->GetPosition(pos) ;
457 p->SetFiducialArea(GetFiducialArea(pos)) ;
458\r
459 fhRecAll->Fill( p->Pt() ) ; //All recontructed particles\r
460 Int_t iFidArea = p->GetFiducialArea(); \r
461 if(iFidArea>0){\r
462 fhRecAllArea1->Fill(p->Pt() ) ; \r
463 if(iFidArea>1){\r
464 fhRecAllArea2->Fill(p->Pt() ) ; \r
465 if(iFidArea>2){\r
466 fhRecAllArea3->Fill(p->Pt() ) ; \r
467 }\r
468 }\r
469 }\r
470\r
471 if(fStack){\r
472 TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;\r
473 if(fDebug>2)\r
474 printf("Pdgcode = %d\n",prim->GetPdgCode());\r
475\r
476 if(prim->GetPdgCode()!=22){ //not photon\r
477//<--DP p->SetPhoton(kFALSE);\r
478 fhRecOther->Fill(p->Pt()); //not photons real spectra\r
479 for(Int_t iPID=0; iPID<4; iPID++){\r
480 if(p->IsPIDOK(iPID,22))\r
481 fhRecOtherPID[iPID]->Fill(p->Pt());\r
482 }\r
483 }\r
484 else{ //Primary photon (as in MC)\r
485//<--DP p->SetPhoton(kTRUE);\r
486 fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon\r
487 for(Int_t iPID=0; iPID<4; iPID++){\r
488 if(p->IsPIDOK(iPID,22))\r
489 fhRecPhotonPID[iPID]->Fill(p->Pt());\r
490 }\r
491 Int_t pi0i=prim->GetFirstMother();\r
492 Int_t grandpaPDG=-1 ;\r
493 TParticle * pi0p = 0;\r
494 if(pi0i>=0){\r
495 pi0p=fStack->Particle(pi0i);\r
496 grandpaPDG=pi0p->GetPdgCode() ;\r
497 }\r
498 switch(grandpaPDG){\r
499 case 111: //Pi0 decay\r
500 //Primary decay photon (as in MC)\r
501 fhRecPhotPi0->Fill(p->Pt());\r
502 break ;\r
503 case 11:\r
504 case -11: //electron/positron conversion\r
505//<--DP p->SetConverted(1);\r
506 fhRecPhotConv->Fill(p->Pt()); //Reconstructed with photon from conversion primary\r
507 fhConversionRadius->Fill(prim->R());\r
508 break ;\r
509 case -2212:\r
510 case -2112: //antineutron & antiproton conversion\r
511 fhRecPhotHadron->Fill(p->Pt()); //Reconstructed with photon from antibaryon annihilation\r
512 fhInteractionRadius->Fill(prim->R());\r
513 break ;\r
514 \r
515 case 221: //eta decay\r
516 fhRecPhotEta->Fill(p->Pt());\r
517 break ; \r
518 \r
519 case 223: //omega meson decay\r
520 fhRecPhotOmega->Fill(p->Pt());\r
521 break ;\r
522 \r
523 case 331: //eta' decay\r
524 fhRecPhotEtapr->Fill(p->Pt());\r
525 break ;\r
526 \r
527 case -1: //direct photon or no primary\r
528 fhRecPhotDirect->Fill(p->Pt());\r
529 break ;\r
530 \r
531 default: \r
532 fhRecPhotOther->Fill(p->Pt());\r
533 break ;\r
534 } \r
535\r
536 //Now classify pi0 decay photon\r
537 if(grandpaPDG==111){\r
538//<--DP p->Pi0Decay(kTRUE); //Mark this photon as primary decayed\r
539//<--DP p->Pi0Id(pi0i); //remember id of the parent pi0\r
540\r
541 //Now check if second (partner) photon from pi0 decay hits PHOS or not\r
542 //i.e. both photons can be tagged or it's the systematic error\r
543//<--DP p->SetPartnerPt(0.); \r
544 Int_t indexdecay1,indexdecay2;\r
545\r
546 indexdecay1=pi0p->GetFirstDaughter();\r
547 indexdecay2=pi0p->GetLastDaughter();\r
548 Int_t indexdecay=-1;\r
549 if(fDebug>2)\r
550 printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());\r
551 \r
552 if(indexdecay1!=caloCluster->GetLabel()) \r
553 indexdecay=indexdecay1;\r
554 if(indexdecay2!=caloCluster->GetLabel()) \r
555 indexdecay=indexdecay2;\r
556 if(indexdecay==-1){\r
557 if(fDebug>2){\r
558 printf("Probably the other photon is not in the stack!\n");\r
559 printf("Number of daughters: %d\n",pi0p->GetNDaughters());\r
560 }\r
561 } \r
562 else{\r
563 TParticle *partner = fStack->Particle(indexdecay);\r
564 Int_t modulenum;\r
565 Double_t xtmp,ztmp;\r
566//<--DP p->SetPartnerPt(partner->Pt());\r
567 if(partner->GetPdgCode()==22){ \r
568 Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason\r
569 if(partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS\r
570 if(fDebug>2)\r
571 printf("P_Reg, E=%f\n",partner->Energy());\r
572 fhPartnerMissedEmin->Fill(partner->Pt()); //Spectrum of missed partners\r
573 fhDecWMissedPartnerEmin->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
574 isPartnerLost=kTRUE;\r
575 }\r
576 if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector\r
577 if(fDebug>2)\r
578 printf("P_Conv, daughters=%d\n",partner->GetNDaughters());\r
579//<--DP p->SetConvertedPartner(1);\r
580 fhPartnerMissedConv->Fill(partner->Pt());\r
581 fhDecWMissedPartnerConv->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
582 isPartnerLost=kTRUE;\r
583 }\r
584 if(!fgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp)){ //this photon cannot hit PHOS\r
585 if(fDebug>2)\r
586 printf("P_Geo\n");\r
587 fhPartnerMissedGeo->Fill(partner->Pt());\r
588 fhDecWMissedPartnerGeom0->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
589 if(iFidArea>0){\r
590 fhDecWMissedPartnerGeom1->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
591 if(iFidArea>1){\r
592 fhDecWMissedPartnerGeom2->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
593 if(iFidArea>2){\r
594 fhDecWMissedPartnerGeom3->Fill(p->Pt()) ; //Spectrum of tagged with missed partner\r
595 }\r
596 }\r
597 }\r
598 isPartnerLost=kTRUE;\r
599 }\r
600 if(!isPartnerLost){\r
601// p->SetMCTagged(1); //set this photon as primary tagged\r
602 fhDecWMCPartner->Fill(p->Pt());\r
603 fhPartnerMCReg->Fill(partner->Pt());\r
604 if(fDebug>2){\r
605 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", caloCluster->E(),partner->Pt(),partner->Energy(),partner->Px(),partner->Py(),partner->Pz(),partner->GetNDaughters());\r
606 printf(", module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);\r
607 }\r
608 }\r
609 else{\r
610 fhDecWMissedPartnerAll->Fill(p->Pt());\r
611 }\r
612 }//Partner - photon\r
613 else{//partner not photon\r
614 fhDecWMissedPartnerNotPhoton->Fill(p->Pt()); \r
615 }\r
616 }\r
617 }\r
618 }\r
619 }\r
620 } //PHOS clusters\r
621 \r
622 if(fDebug>1) \r
623 printf("number of clusters: %d\n",inList);\r
624\r
625 //Invariant Mass analysis\r
626 for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList-1 ; phosPhoton1++) {\r
627 AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1)); \r
628 for(Int_t phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < inList ; phosPhoton2++) {\r
629 AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));\r
630\r
631 Double_t invMass = p1->GetPairMass(p2);\r
632 fhInvMassReal->Fill(invMass,p1->Pt());\r
633 fhInvMassReal->Fill(invMass,p2->Pt());\r
634 if(fDebug>2)\r
635 printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);\r
636\r
637 Bool_t makePi01=IsInPi0Band(invMass,p1->Pt());\r
638 Bool_t makePi02=IsInPi0Band(invMass,p2->Pt());\r
639\r
640\r
641 if(makePi01 && p1->IsTagged()){//Multiple tagging\r
642 fhTaggedMult->Fill(p1->Pt());\r
643 } \r
644 if(makePi01 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times\r
645 fhTaggedAll->Fill(p1->Pt());\r
646 Int_t iFidArea = p1->GetFiducialArea(); \r
647 if(iFidArea>0){\r
648 fhTaggedArea1->Fill(p1->Pt() ) ; \r
649 if(iFidArea>1){\r
650 fhTaggedArea2->Fill(p1->Pt() ) ; \r
651 if(iFidArea>2){\r
652 fhTaggedArea3->Fill(p1->Pt() ) ; \r
653 }\r
654 }\r
655 }\r
656\r
657 for(Int_t iPID=0; iPID<4; iPID++){\r
658 if(p1->IsPIDOK(iPID,22))\r
659 fhTaggedPID[iPID]->Fill(p1->Pt());\r
660 }\r
661\r
662 p1->SetTagged(kTRUE) ;\r
663 } \r
664 if(makePi02 && p2->IsTagged()){//Multiple tagging\r
665 fhTaggedMult->Fill(p2->Pt());\r
666 } \r
667 if(makePi02 && !p2->IsTagged()){//How should be account for multiply tagged photons?\r
668 fhTaggedAll->Fill(p2->Pt());\r
669 p2->SetTagged(kTRUE) ;\r
670 }\r
671 \r
672 //Now get use MC information\r
673 //First chesk if this is true pi0 pair\r
674 if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0\r
675// p1->SetTrueTagged(1);\r
676// p2->SetTrueTagged(1);\r
677 if(makePi01)//Correctly tagged photons\r
678 fhTaggedMCTrue->Fill(p1->Pt());\r
679 else{ //Decay pair missed tagging \r
680 fhMCMissedTagging->Fill(p1->Pt());\r
681 fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;\r
682 //Clussify why missed tagging (todo)\r
683 //Converted\r
684 //Partner not a photon\r
685 //Tagged not a photon\r
686 //Just wrong inv.mass \r
687 } \r
688 if(makePi02)\r
689 fhTaggedMCTrue->Fill(p2->Pt());\r
690 else{ \r
691 fhMCMissedTagging->Fill(p2->Pt());\r
692 fhMCMissedTaggingMass->Fill(invMass,p2->Pt()) ;\r
693 //Clussify why missed tagging (todo)\r
694 //Converted\r
695 //Partner not a photon\r
696 //Tagged not a photon\r
697 //Just wrong inv.mass \r
698 } \r
699 }\r
700 else{//Fake tagged - not from the same pi0\r
701 if(makePi01)//Fake pair\r
702 fhMCFakeTagged->Fill(p1->Pt());\r
703 if(makePi02)\r
704 fhMCFakeTagged->Fill(p2->Pt());\r
705 }\r
706 }\r
707 }\r
708\r
709 //Fill Mixed InvMass distributions:\r
710 TIter nextEv(fEventList) ;\r
711 while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){\r
712 Int_t nPhotons2 = event2->GetEntriesFast() ;\r
713 for(Int_t i=0; i < inList ; i++){\r
714 AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;\r
715 for(Int_t j=0; j < nPhotons2 ; j++){\r
716 AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;\r
717 Double_t invMass = p1->GetPairMass(p2);\r
718 fhInvMassMixed->Fill(invMass,p1->Pt());\r
719 fhInvMassMixed->Fill(invMass,p2->Pt());\r
720 }\r
721 }\r
722 }\r
723\r
724/* if(inList==1){ //We have only one photon in PHOS, but still there can be missing partner for it\r
725 AliAODPWG4Particle * p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(0));\r
726// if(p->IsMCTagged()){\r
727 //should have pair but partner was not registered in PHOS\r
728// double pairpt=p->GetPairPt();\r
729// if(fDebug>2)\r
730// printf("pair not found for phot: E=%f, pair Pt=%f\n",p->Energy(),pairpt);\r
731// fhMCPartnerMissed->Fill(pairpt);\r
732// fhstra->Fill(p->Pt());\r
733// }\r
734 }\r
735*/\r
736\r
737 //Remove old events\r
738 fEventList->AddFirst(fCaloPhotonsArr);\r
739 if(fEventList->GetSize() > 10){\r
740 TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());\r
741 fEventList->Remove(tmp);\r
742 delete tmp;\r
743 }\r
744\r
745 PostData(1, fOutputList);\r
746}\r
747\r
748\r
749//______________________________________________________________________________\r
750void AliAnalysisTaskTaggedPhotons::Init()\r
751{\r
752 // Intialisation of parameters\r
753 AliInfo("Doing initialisation") ; \r
754 SetPhotonId(0.9) ; \r
755 SetMinEnergyCut(0.4);\r
756 SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);\r
757 SetPi0SigmaParameters(0.004508,0.005497,0.00000006);\r
758}\r
759\r
760//______________________________________________________________________________\r
761void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)\r
762{\r
763 // Processing when the event loop is ended\r
764\r
765 //Write everything to the file\r
766 char outname[55];\r
767 if(fPHOS)\r
768 sprintf(outname,"Tagging_PHOS.root") ;\r
769 else \r
770 sprintf(outname,"Tagging_EMCAL.root") ;\r
771 TFile *outfile = new TFile (outname,"recreate");\r
772\r
773fhRecAll->Write();\r
774fhRecAllArea1->Write();\r
775fhRecAllArea2->Write();\r
776fhRecAllArea3->Write();\r
777fhRecPhoton->Write();\r
778fhRecOther->Write();\r
779fhRecPhotonPID[0]->Write();\r
780fhRecPhotonPID[1]->Write();\r
781fhRecPhotonPID[2]->Write();\r
782fhRecPhotonPID[3]->Write();\r
783fhRecOtherPID[0]->Write();\r
784fhRecOtherPID[1]->Write();\r
785fhRecOtherPID[2]->Write();\r
786fhRecOtherPID[3]->Write();\r
787fhRecPhotPi0->Write();\r
788fhRecPhotEta->Write();\r
789fhRecPhotOmega->Write();\r
790fhRecPhotEtapr->Write();\r
791fhRecPhotConv->Write();\r
792fhRecPhotHadron->Write();\r
793fhRecPhotDirect->Write();\r
794fhRecPhotOther->Write();\r
795fhDecWMCPartner->Write();\r
796fhDecWMissedPartnerNotPhoton->Write();\r
797fhDecWMissedPartnerAll->Write();\r
798fhDecWMissedPartnerEmin->Write();\r
799fhDecWMissedPartnerConv->Write();\r
800fhDecWMissedPartnerGeom0->Write();\r
801fhDecWMissedPartnerGeom1->Write();\r
802fhDecWMissedPartnerGeom2->Write();\r
803fhDecWMissedPartnerGeom3->Write();\r
804fhPartnerMCReg->Write();\r
805fhPartnerMissedEmin->Write();\r
806fhPartnerMissedConv->Write();\r
807fhPartnerMissedGeo->Write();\r
808fhTaggedAll->Write();\r
809fhTaggedArea1->Write();\r
810fhTaggedArea2->Write();\r
811fhTaggedArea3->Write();\r
812\r
813fhTaggedPID[0]->Write();\r
814fhTaggedPID[1]->Write();\r
815fhTaggedPID[2]->Write();\r
816fhTaggedPID[3]->Write();\r
817fhTaggedMult->Write();\r
818fhTaggedMCTrue->Write();\r
819fhMCMissedTagging->Write();\r
820fhMCFakeTagged->Write();\r
821fhInvMassReal->Write();\r
822fhInvMassMixed->Write();\r
823fhMCMissedTaggingMass->Write();\r
824fhConversionRadius->Write();\r
825fhInteractionRadius->Write();\r
826fhEvents->Write();\r
827\r
828/*\r
829 fhAllPhotons->Write() ;\r
830 fhNotPhotons->Write() ;\r
831 fhAllPhotonsPrimary->Write() ;\r
832 fhNotPhotonsPrimary->Write() ;\r
833 fhfakeNotPhotons->Write() ;\r
834\r
835 fhTaggedPhotons->Write();\r
836 fhfakeTaggedPhotons->Write();\r
837 fhDecayNotTaggedPhotons->Write();\r
838 fhstrangeNotTaggedPhotons->Write();\r
839 fhstrangeNotTaggedPhotonsPair->Write();\r
840 fhstrangeNotTaggedPhotonsRegCut->Write();\r
841 fhstrangeNotTaggedPhotonsPairRegCut->Write();\r
842\r
843 fhPi0DecayPhotonsPrimary->Write();\r
844 fhEtaDecayPhotonsPrimary->Write();\r
845 fhOmegaDecayPhotonsPrimary->Write();\r
846 fhEtaSDecayPhotonsPrimary->Write();\r
847 fhOtherDecayPhotonsPrimary->Write();\r
848 fhDecayPhotonsPrimary->Write();\r
849 fhConvertedPhotonsPrimary->Write();\r
850 fhConvertedPhotonsPrimaryHadronsDecays->Write();\r
851 fhCoordsConvertion->Write();\r
852 fhCoordsConvertion2->Write();\r
853\r
854 fhPHOSInvariantMassReal->Write();\r
855 fhPHOSInvariantMassMixed->Write();\r
856\r
857 fhPHOSPi0->Write();\r
858\r
859 fhPi0DecayPhotonsGeomfake->Write();\r
860 fhPi0DecayPhotonsTaggedPrimary->Write();\r
861 fhPi0DecayPhotonsTaggedPrimaryPair->Write();\r
862 fhPi0DecayPhotonsTaggedPrimaryRegCut->Write();\r
863 fhPi0DecayPhotonsTaggedPrimaryPairRegCut->Write();\r
864\r
865 fhPi0DecayPhotonsBigDecay->Write();\r
866 fhPi0DecayPhotonsPConv->Write();\r
867 fhPi0DecayPhotonsPGeo->Write();\r
868 fhPi0DecayPhotonsPReg->Write();\r
869\r
870 fhfakeTaggedPhotonsConv->Write();\r
871 fhfakeTaggedPhotonsPID->Write();\r
872\r
873 fhTrackRefCoords->Write();\r
874 fhEvents->Write();\r
875*/\r
876\r
877outfile->Close();\r
878\r
879}\r
880//______________________________________________________________________________\r
881Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)\r
882{\r
883 //Parameterization of the pi0 peak region\r
884 Double_t mpi0_mean = fPi0Mean_p0 + fPi0Mean_p1 * pt + fPi0Mean_p2 * pt*pt + fPi0Mean_p3 * pt*pt*pt;\r
885 Double_t mpi0_sigma = TMath::Sqrt(fPi0Sigma_p0 * fPi0Sigma_p0 / pt + fPi0Sigma_p1 * fPi0Sigma_p1 + fPi0Sigma_p2 * fPi0Sigma_p2 / pt / pt);\r
886 \r
887 return (m>mpi0_mean-2*mpi0_sigma && m>mpi0_mean+2*mpi0_sigma) ;\r
888}\r
889//______________________________________________________________________________\r
890Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(AliAODPWG4Particle *p1, AliAODPWG4Particle *p2){\r
891 //Looks through parents and finds if there was commont pi0 among ancestors\r
892\r
893 if(!fStack)\r
894 return kFALSE ; //can not say anything\r
895\r
896 Int_t prim1 = p1->GetLabel();\r
897 while(prim1!=-1){ \r
898 Int_t prim2 = p2->GetLabel();\r
899 while(prim2!=-1){ \r
900 if(prim1==prim2){\r
901 if(fStack->Particle(prim1)->GetPdgCode()==111)\r
902 return kTRUE ;\r
903 else\r
904 return kFALSE ;\r
905 }\r
906 prim2=fStack->Particle(prim2)->GetFirstMother() ;\r
907 }\r
908 prim1=fStack->Particle(prim1)->GetFirstMother() ;\r
909 }\r
910 return kFALSE ;\r
911}\r
912//______________________________________________________________________________\r
913Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos){\r
914 //calculates in which kind of fiducial area photon hit\r
915 if(fPHOS){\r
916 Double_t phi=TMath::ATan2(pos[1],pos[0]) ;\r
917 Double_t z=pos[2] ;
918 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;\r
919 while(phi<0.)phi+=TMath::TwoPi() ;\r
920 //From active PHOS areat remove bands in 10 cm\r
921 const Double_t kphi=TMath::ATan(10./460.) ; //angular band width\r
922 Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;\r
923 Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;\r
924 Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);\r
925 Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);\r
926 return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); \r
927 }\r
928 else{//EMCAL\r
929 //For the moment whole EMCAL is considered as area 1\r
930 return 1 ;\r
931 }\r
932}\r