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