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