]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaOmegaToPi0Gamma.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
correct binning of hNCellsMod histogram, set remaining TH3 histogram to be filled...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaOmegaToPi0Gamma.cxx
CommitLineData
21a4b1c0 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/* $Id: $ */
16//_________________________________________________________________________
17// class to extract omega(782)->pi0+gamma->3gamma
18// Mar. 22, 2011: Additional method, espeically for EMCAL. A high E cluster is assumpted as pi0 (two photons are overlapped) without unfolding
19//
20//-- Author: Renzhuo Wan (IOPP-Wuhan, China)
21//_________________________________________________________________________
22
23// --- ROOT system
24class TROOT;
25
26// --- AliRoot system
27//class AliVEvent;
28// --- ROOT system ---
29#include "TH2F.h"
30#include "TLorentzVector.h"
31#include "TParticle.h"
32#include "TCanvas.h"
33#include "TFile.h"
34//---- AliRoot system ----
35#include "AliAnaOmegaToPi0Gamma.h"
36#include "AliCaloTrackReader.h"
37#include "AliCaloPID.h"
38#include "AliStack.h"
39#include "AliVEvent.h"
40#include "AliAODEvent.h"
41#include "AliAODMCParticle.h"
42ClassImp(AliAnaOmegaToPi0Gamma)
43
44//______________________________________________________________________________
45AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma() : AliAnaPartCorrBaseClass(),
46fInputAODPi0(0), fInputAODGammaName(""),
47fEventsList(0x0),fNVtxZBin(0), fNCentBin(0), fNRpBin(0), fNBadChDistBin(0), fNpid(0),
48fNmaxMixEv(0), fVtxZCut(0), fCent(0), fRp(0),
49fPi0Mass(0),fPi0MassWindow(0),fPi0OverOmegaPtCut(0),
50fGammaOverOmegaPtCut(0), fEOverlapCluster(0),
51fhEtalon(0),
52fRealOmega0(0), fMixAOmega0(0),
53fMixBOmega0(0), fMixCOmega0(0),
54fRealOmega1(0), fMixAOmega1(0),
55fMixBOmega1(0), fMixCOmega1(0),
56fRealOmega2(0), fMixAOmega2(0),
57fMixBOmega2(0), fMixCOmega2(0),
58fhFakeOmega(0),
59fhOmegaPriPt(0)
60{
61 //Default Ctor
62 InitParameters();
63}
64/*
65//______________________________________________________________________________
66AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma(const AliAnaOmegaToPi0Gamma & ex) : AliAnaPartCorrBaseClass(ex),
67fInputAODPi0(new TClonesArray (*ex.fInputAODPi0)),
68fInputAODGammaName(ex.fInputAODGammaName),
69fEventsList(0x0),
70fNVtxZBin(ex.fNVtxZBin), fNCentBin(ex.fNCentBin), fNRpBin(ex.fNRpBin),
71fNBadChDistBin(ex.fNBadChDistBin),fNpid(ex.fNpid),
72fNmaxMixEv(ex.fNmaxMixEv),
73fVtxZCut(ex.fVtxZCut), fCent(ex.fCent), fRp(ex.fRp),
74fPi0Mass(ex.fPi0Mass),
75fPi0MassWindow(ex.fPi0MassWindow),
76fPi0OverOmegaPtCut(ex.fPi0OverOmegaPtCut),
77fGammaOverOmegaPtCut(ex.fGammaOverOmegaPtCut),
78fhEtalon(ex.fhEtalon),
79fRealOmega0(ex.fRealOmega0), fMixAOmega0(ex.fMixAOmega0),
80fMixBOmega0(ex.fMixBOmega0), fMixCOmega0(ex.fMixCOmega0),
81fRealOmega1(ex.fRealOmega1), fMixAOmega1(ex.fMixAOmega1),
82fMixBOmega1(ex.fMixBOmega1), fMixCOmega1(ex.fMixCOmega1),
83fRealOmega2(ex.fRealOmega2), fMixAOmega2(ex.fMixAOmega2),
84fMixBOmega2(ex.fMixBOmega2), fMixCOmega2(ex.fMixCOmega2),
85fhOmegaPriPt(ex.fhOmegaPriPt)
86{
87 // cpy ctor
88 //Do not need it
89}
90*/
91/*
92//______________________________________________________________________________
93AliAnaOmegaToPi0Gamma & AliAnaOmegaToPi0Gamma::operator = (const AliAnaOmegaToPi0Gamma & ex)
94{
95 // assignment operator
96
97 if(this == &ex)return *this;
98 ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
99 fInputAODPi0 = new TClonesArray(*ex.fInputAODPi0);
100 fInputAODGammaName = ex.fInputAODGammaName;
101 fEventsList = ex.fEventsList;
102
103 fNVtxZBin=ex.fNVtxZBin;
104 fNCentBin=ex.fNCentBin;
105 fNRpBin=ex.fNRpBin;
106 fNBadChDistBin=ex.fNBadChDistBin;
107 fNpid=ex.fNpid;
108 fNmaxMixEv =ex.fNmaxMixEv;
109
110 fVtxZCut=ex.fVtxZCut;
111 fCent=ex.fCent;
112 fRp=ex.fRp;
113
114 fPi0Mass=ex.fPi0Mass;
115 fPi0MassWindow=ex.fPi0MassWindow;
116 fPi0OverOmegaPtCut=ex.fPi0OverOmegaPtCut;
117 fGammaOverOmegaPtCut=ex.fGammaOverOmegaPtCut;
118
119 fhEtalon=ex.fhEtalon;
120 fRealOmega0=ex.fRealOmega0;
121 fMixAOmega0=ex.fMixAOmega0;
122 fMixBOmega0=ex.fMixBOmega0;
123 fMixCOmega0=ex.fMixCOmega0;
124 fRealOmega1=ex.fRealOmega1;
125 fMixAOmega1=ex.fMixAOmega1;
126 fMixBOmega1=ex.fMixBOmega1;
127 fMixCOmega1=ex.fMixCOmega1;
128 fRealOmega2=ex.fRealOmega2;
129 fMixAOmega2=ex.fMixAOmega2;
130 fMixBOmega2=ex.fMixBOmega2;
131 fMixCOmega2=ex.fMixCOmega2;
132 fhOmegaPriPt=ex.fhOmegaPriPt;
133 return *this;
134
135}
136*/
137
138//______________________________________________________________________________
139AliAnaOmegaToPi0Gamma::~AliAnaOmegaToPi0Gamma() {
140
141 //dtor
142// Done by the maker
143// if(fInputAODPi0){
144// fInputAODPi0->Clear();
145// delete fInputAODPi0;
146// }
147
148 if(fEventsList){
149 for(Int_t i=0;i<fNVtxZBin;i++){
150 for(Int_t j=0;j<fNCentBin;j++){
151 for(Int_t k=0;k<fNRpBin;k++){
152 fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->Clear();
153 delete fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k];
154 }
155 }
156 }
157 }
158 delete [] fEventsList;
159 fEventsList=0;
160
161 delete [] fVtxZCut;
162 delete [] fCent;
163 delete [] fRp;
164
165}
166
167//______________________________________________________________________________
168void AliAnaOmegaToPi0Gamma::InitParameters()
169{
170//Init parameters when first called the analysis
171//Set default parameters
172 fInputAODGammaName="PhotonsDetector";
173 fNVtxZBin=1;
174 fNCentBin=1;
175 fNRpBin=1;
176 fNBadChDistBin=3;
177 fNpid=1;
178 fNmaxMixEv=8;
179
180 fPi0Mass=0.1348;
181 fPi0MassWindow=0.015;
182 fPi0OverOmegaPtCut=0.8;
183 fGammaOverOmegaPtCut=0.2;
184
185 fEOverlapCluster=6;
186}
187
188
189//______________________________________________________________________________
190TList * AliAnaOmegaToPi0Gamma::GetCreateOutputObjects()
191{
192 //
193 fVtxZCut = new Double_t [fNVtxZBin];
194 for(Int_t i=0;i<fNVtxZBin;i++) fVtxZCut[i]=10*(i+1);
195
196 fCent=new Double_t[fNCentBin];
197 for(int i = 0;i<fNCentBin;i++)fCent[i]=0;
198
199 fRp=new Double_t[fNRpBin];
200 for(int i = 0;i<fNRpBin;i++)fRp[i]=0;
201 //
202 Int_t nptbins = GetHistoPtBins();
203 Float_t ptmax = GetHistoPtMax();
204 Float_t ptmin = GetHistoPtMin();
205
206 Int_t nmassbins = GetHistoMassBins();
207 Float_t massmin = GetHistoMassMin();
208 Float_t massmax = GetHistoMassMax();
209
210 fhEtalon = new TH2F("hEtalon","Histo with binning parameters", nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
211 fhEtalon->SetXTitle("P_{T} (GeV)") ;
212 fhEtalon->SetYTitle("m_{inv} (GeV)") ;
213
214 // store them in fOutputContainer
215 fEventsList = new TList*[fNVtxZBin*fNCentBin*fNRpBin];
216 for(Int_t i=0;i<fNVtxZBin;i++){
217 for(Int_t j=0;j<fNCentBin;j++){
218 for(Int_t k=0;k<fNRpBin;k++){
219 fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k] = new TList();
220 fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->SetOwner(kFALSE);
221 }
222 }
223 }
224
225 TList * outputContainer = new TList() ;
226 outputContainer->SetName(GetName());
227 const Int_t buffersize = 255;
228 char key[buffersize] ;
229 char title[buffersize] ;
230 const char * detector= fInputAODGammaName.Data();
231 Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
232
233 fRealOmega0 =new TH2F*[ndim];
234 fMixAOmega0 =new TH2F*[ndim];
235 fMixBOmega0 =new TH2F*[ndim];
236 fMixCOmega0 =new TH2F*[ndim];
237
238 fRealOmega1 =new TH2F*[ndim];
239 fMixAOmega1 =new TH2F*[ndim];
240 fMixBOmega1 =new TH2F*[ndim];
241 fMixCOmega1 =new TH2F*[ndim];
242
243 fRealOmega2 =new TH2F*[ndim];
244 fMixAOmega2 =new TH2F*[ndim];
245 fMixBOmega2 =new TH2F*[ndim];
246 fMixCOmega2 =new TH2F*[ndim];
247
248 fhFakeOmega = new TH2F*[fNCentBin];
249
250 for(Int_t i=0;i<fNVtxZBin;i++){
251 for(Int_t j=0;j<fNCentBin;j++){
252 for(Int_t k=0;k<fNRpBin;k++){ //at event level
253 Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
254 for(Int_t ipid=0;ipid<fNpid;ipid++){
255 for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle level
256
257 Int_t index=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
258
259 snprintf(key,buffersize,"RealToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
260 snprintf(title,buffersize, "%s Real Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
261 fRealOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
262 fRealOmega0[index]->SetName(key) ;
263 fRealOmega0[index]->SetTitle(title);
264 outputContainer->Add(fRealOmega0[index]);
265
266 snprintf(key,buffersize,"MixAToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
267 snprintf(title,buffersize, "%s MixA Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
268 fMixAOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
269 fMixAOmega0[index]->SetName(key) ;
270 fMixAOmega0[index]->SetTitle(title);
271 outputContainer->Add(fMixAOmega0[index]);
272
273 snprintf(key,buffersize,"MixBToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
274 snprintf(title,buffersize, "%s MixB Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
275 fMixBOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
276 fMixBOmega0[index]->SetName(key) ;
277 fMixBOmega0[index]->SetTitle(title);
278 outputContainer->Add(fMixBOmega0[index]);
279
280 snprintf(key,buffersize,"MixCToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
281 snprintf(title,buffersize, "%s MixC Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
282 fMixCOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
283 fMixCOmega0[index]->SetName(key) ;
284 fMixCOmega0[index]->SetTitle(title);
285 outputContainer->Add(fMixCOmega0[index]);
286
287 snprintf(key,buffersize,"RealToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
288 snprintf(title,buffersize, "%s Real Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
289 fRealOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
290 fRealOmega1[index]->SetName(key) ;
291 fRealOmega1[index]->SetTitle(title);
292 outputContainer->Add(fRealOmega1[index]);
293
294 snprintf(key,buffersize,"MixAToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
295 snprintf(title,buffersize, "%s MixA Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
296 fMixAOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
297 fMixAOmega1[index]->SetName(key) ;
298 fMixAOmega1[index]->SetTitle(title);
299 outputContainer->Add(fMixAOmega1[index]);
300
301 snprintf(key,buffersize,"MixBToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
302 snprintf(title,buffersize, "%s MixB Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
303 fMixBOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
304 fMixBOmega1[index]->SetName(key) ;
305 fMixBOmega1[index]->SetTitle(title);
306 outputContainer->Add(fMixBOmega1[index]);
307
308 snprintf(key,buffersize,"MixCToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
309 snprintf(title,buffersize, "%s MixC Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
310 fMixCOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
311 fMixCOmega1[index]->SetName(key) ;
312 fMixCOmega1[index]->SetTitle(title);
313 outputContainer->Add(fMixCOmega1[index]);
314
315 snprintf(key,buffersize,"RealToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
316 snprintf(title,buffersize, "%s Real Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
317 fRealOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
318 fRealOmega2[index]->SetName(key) ;
319 fRealOmega2[index]->SetTitle(title);
320 outputContainer->Add(fRealOmega2[index]);
321
322 snprintf(key,buffersize,"MixAToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
323 snprintf(title,buffersize, "%s MixA Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
324 fMixAOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
325 fMixAOmega2[index]->SetName(key) ;
326 fMixAOmega2[index]->SetTitle(title);
327 outputContainer->Add(fMixAOmega2[index]);
328
329 snprintf(key,buffersize,"MixBToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
330 snprintf(title,buffersize, "%s MixB Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
331 fMixBOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
332 fMixBOmega2[index]->SetName(key) ;
333 fMixBOmega2[index]->SetTitle(title);
334 outputContainer->Add(fMixBOmega2[index]);
335
336 snprintf(key,buffersize,"MixCToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
337 snprintf(title,buffersize, "%s MixC Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
338 fMixCOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
339 fMixCOmega2[index]->SetName(key) ;
340 fMixCOmega2[index]->SetTitle(title);
341 outputContainer->Add(fMixCOmega2[index]);
342 }
343 }
344 }
345 }
346 }
347
348 for(Int_t i=0;i<fNCentBin;i++){
349 snprintf(key,buffersize, "fhFakeOmega%d",i);
350 snprintf(title,buffersize,"FakePi0(high pt cluster)+Gamma with Centrality%d",i);
351 fhFakeOmega[i] = (TH2F*)fhEtalon->Clone(key);
352 fhFakeOmega[i]->SetTitle(title);
353 outputContainer->Add(fhFakeOmega[i]);
354 }
355 if(IsDataMC()){
356 snprintf(key,buffersize, "%sOmegaPri",detector);
357 snprintf(title,buffersize,"primary #omega in %s",detector);
358 fhOmegaPriPt=new TH1F(key, title,nptbins,ptmin,ptmax);
359 fhOmegaPriPt->GetXaxis()->SetTitle("P_{T}");
360 fhOmegaPriPt->GetYaxis()->SetTitle("dN/P_{T}");
361 outputContainer->Add(fhOmegaPriPt);
362 }
363
364 delete fhEtalon;
365 return outputContainer;
366}
367
368//______________________________________________________________________________
369void AliAnaOmegaToPi0Gamma::Print(const Option_t * /*opt*/) const
370{
371 //Print some relevant parameters set in the analysis
372 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
373 AliAnaPartCorrBaseClass::Print(" ");
374 printf("Omega->pi0+gamma->3gamma\n");
375 printf("Cuts at event level: \n");
376 printf("Bins of vertex Z: %d \n", fNVtxZBin);
377 printf("Bins of centrality: %d \n",fNCentBin);
378 printf("Bins of Reaction plane: %d\n",fNRpBin);
379 printf("Cuts at AOD particle level:\n");
380 printf("Number of PID: %d \n", fNpid);
381 printf("Number of DistToBadChannel cuts: %d\n", fNBadChDistBin);
382 printf("number of events buffer to be mixed: %d\n",fNmaxMixEv);
383}
384
385//______________________________________________________________________________
386void AliAnaOmegaToPi0Gamma::MakeAnalysisFillHistograms()
387{
388 //fill the MC AOD if needed first
389 //-----------
390 //need to be further implemented
391 AliStack * stack = 0x0;
392 // TParticle * primary = 0x0;
393 TClonesArray * mcparticles0 = 0x0;
394 //TClonesArray * mcparticles1 = 0x0;
395 AliAODMCParticle * aodprimary = 0x0;
396 Int_t pdg=0;
397 Double_t pt=0;
398 Double_t eta=0;
399
400 if(IsDataMC()){
401 if(GetReader()->ReadStack()){
402 stack = GetMCStack() ;
403 if(!stack){
404 printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - There is no stack!\n");
405 }
406 else{
407 for(Int_t i=0 ; i<stack->GetNtrack(); i++){
408 TParticle * prim = stack->Particle(i) ;
409 pdg = prim->GetPdgCode() ;
410 eta=prim->Eta();
411 pt=prim->Pt();
412 if(TMath::Abs(eta)<0.5) {
413 if(pdg==223) fhOmegaPriPt->Fill(pt);
414 }
415 }
416 }
417 }
418 else if(GetReader()->ReadAODMCParticles()){
419 //Get the list of MC particles
420 mcparticles0 = GetReader()->GetAODMCParticles(0);
421 if(!mcparticles0 ) {
422 if(GetDebug() > 0) printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
423 }
424 // if(GetReader()->GetSecondInputAODTree()){
425 // mcparticles1 = GetReader()->GetAODMCParticles(1);
426 // if(!mcparticles1 && GetDebug() > 0) {
427 // printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
428 // }
429 // }
430 else{
431 for(Int_t i=0;i<mcparticles0->GetEntries();i++){
432 aodprimary =(AliAODMCParticle*)mcparticles0->At(i);
433 pdg = aodprimary->GetPdgCode() ;
434 eta=aodprimary->Eta();
435 pt=aodprimary->Pt();
436 if(TMath::Abs(eta)<0.5) {
437 if(pdg==223) fhOmegaPriPt->Fill(pt);
438 }
439
440 }
441 }//mcparticles0 exists
442 }//AOD MC Particles
443 }// is data and MC
444
445
446 //process event from AOD brach
447 //extract pi0, eta and omega analysis
448 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
449 if(IsBadRun(iRun)) return ;
450
451 //vertex z
452 Double_t vert[]={0,0,0} ;
453 GetVertex(vert);
454 Int_t curEventBin =0;
455
456 Int_t ivtxzbin=(Int_t)TMath::Abs(vert[2])/10;
457 if(ivtxzbin>=fNVtxZBin)return;
458
459 //centrality
460 Int_t currentCentrality = GetEventCentrality();
461 if(currentCentrality == -1) return;
462 Int_t optCent = GetReader()->GetCentralityOpt();
463 Int_t icentbin=currentCentrality/(optCent/fNCentBin) ; //GetEventCentrality();
464
465 printf("-------------- %d %d %d ",currentCentrality, optCent, icentbin);
466 //reaction plane
467 Int_t irpbin=0;
468
469 if(ivtxzbin==-1) return;
470 curEventBin = ivtxzbin*fNCentBin*fNRpBin + icentbin*fNRpBin + irpbin;
471 TClonesArray *aodGamma = (TClonesArray*) GetAODBranch(fInputAODGammaName); //photon array
472 // TClonesArray *aodGamma = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaName); //photon array
473 Int_t nphotons =0;
474 if(aodGamma) nphotons= aodGamma->GetEntries();
475 else return;
476
477 fInputAODPi0 = (TClonesArray*)GetInputAODBranch(); //pi0 array
478 Int_t npi0s = 0;
479 if(fInputAODPi0) npi0s= fInputAODPi0 ->GetEntries();
480 else return;
481
482 if(nphotons<3 || npi0s<1)return; //for pi0, eta and omega->pi0+gamma->3gamma reconstruction
483
484 //reconstruction of omega(782)->pi0+gamma->3gamma
485 //loop for pi0 and photon
486 if(GetDebug() > 0) printf("omega->pi0+gamma->3gamma invariant mass analysis ! This event have %d photons and %d pi0 \n", nphotons, npi0s);
487 for(Int_t i=0;i<npi0s;i++){
488 AliAODPWG4Particle * pi0 = (AliAODPWG4Particle*) (fInputAODPi0->At(i)) ; //pi0
489 TLorentzVector vpi0(pi0->Px(),pi0->Py(),pi0->Pz(),pi0->E());
490 Int_t lab1=pi0->GetCaloLabel(0); // photon1 from pi0 decay
491 Int_t lab2=pi0->GetCaloLabel(1); // photon2 from pi0 decay
492 //for omega->pi0+gamma, it needs at least three photons per event
493 //Get the two decay photons from pi0
494 AliAODPWG4Particle * photon1 =0;
495 AliAODPWG4Particle * photon2 =0;
496 for(Int_t d1=0;d1<nphotons;d1++){
497 for(Int_t d2=0;d2<nphotons;d2++){
498 AliAODPWG4Particle * dp1 = (AliAODPWG4Particle*) (aodGamma->At(d1));
499 AliAODPWG4Particle * dp2 = (AliAODPWG4Particle*) (aodGamma->At(d2));
500 Int_t dlab1=dp1->GetCaloLabel(0);
501 Int_t dlab2=dp2->GetCaloLabel(0);
502 if(dlab1==lab1 && dlab2==lab2){
503 photon1=dp1;
504 photon2=dp2;
505 }
506 else continue;
507 }
508 }
509 //caculate the asy and dist of the two photon from pi0 decay
510 TLorentzVector dph1(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());
511 TLorentzVector dph2(photon2->Px(),photon2->Py(),photon2->Pz(),photon2->E());
512
513 Double_t pi0asy= TMath::Abs(dph1.E()-dph2.E())/(dph1.E()+dph2.E());
514 // Double_t phi1=dph1.Phi();
515 // Double_t phi2=dph2.Phi();
516 // Double_t eta1=dph1.Eta();
517 // Double_t eta2=dph2.Eta();
518 // Double_t pi0dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));
519
520 if(pi0->GetIdentifiedParticleType()==111 && nphotons>2 && npi0s
521 && TMath::Abs(vpi0.M()-fPi0Mass)<fPi0MassWindow) { //pi0 candidates
522
523 //avoid the double counting
524 Int_t * dc1= new Int_t[nphotons];
525 Int_t * dc2= new Int_t[nphotons];
526 Int_t index1=0;
527 Int_t index2=0;
528 for(Int_t k=0;k<i;k++){
529 AliAODPWG4Particle * p3=(AliAODPWG4Particle*)(fInputAODPi0->At(k));
530 Int_t lab4=p3->GetCaloLabel(0);
531 Int_t lab5=p3->GetCaloLabel(1);
532 if(lab1==lab4){ dc1[index1]=lab5; index1++; }
533 if(lab2==lab5){ dc2[index2]=lab4; index2++; }
534 }
535
536
537 //loop the pi0 with third gamma
538 for(Int_t j=0;j<nphotons;j++){
539 AliAODPWG4Particle *photon3 = (AliAODPWG4Particle*) (aodGamma->At(j));
540 TLorentzVector dph3(photon3->Px(),photon3->Py(),photon3->Pz(),photon3->E());
541 Int_t lab3=photon3->GetCaloLabel(0);
542 Double_t pi0gammapt=(vpi0+dph3).Pt();
543 Double_t pi0gammamass=(vpi0+dph3).M();
544 Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
545 Double_t gammaOverOmegaPtRatio= dph3.Pt()/pi0gammapt;
546
547 //pi0, gamma pt cut
548 if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
549 gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
550
551 for(Int_t l=0;l<index1;l++) if(lab3==dc1[l]) lab3=-1;
552 for(Int_t l=0;l<index2;l++) if(lab3==dc2[l]) lab3=-1;
553
554 if(lab3>0 && lab3!=lab1 && lab3!=lab2){
555 for(Int_t ipid=0;ipid<fNpid;ipid++){
556 for(Int_t idist=0;idist<fNBadChDistBin;idist++){
557 Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
558 if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
559 photon2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
560 photon3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
561 photon1->DistToBad()>=idist &&
562 photon2->DistToBad()>=idist &&
563 photon3->DistToBad()>=idist ){
564 //fill the histograms
565 if(GetDebug() > 2) printf("Real: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
566 fRealOmega0[index]->Fill(pi0gammapt,pi0gammamass);
567 if(pi0asy<0.7) fRealOmega1[index]->Fill(pi0gammapt,pi0gammamass);
568 if(pi0asy<0.8) fRealOmega2[index]->Fill(pi0gammapt,pi0gammamass);
569 }
570 }
571 }
572 }
573 }
574 delete []dc1;
575 delete []dc2;
576 if(GetDebug() > 0) printf("MixA: (r1_event1+r2_event1)+r3_event2 \n");
577 //-------------------------
578 //background analysis
579 //three background
580 // --A (r1_event1+r2_event1)+r3_event2
581 Int_t nMixed = fEventsList[curEventBin]->GetSize();
582 for(Int_t im=0;im<nMixed;im++){
583 TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(im));
584 for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
585 AliAODPWG4Particle *mix1ph = (AliAODPWG4Particle*) (ev2->At(mix1));
586 TLorentzVector vmixph(mix1ph->Px(),mix1ph->Py(),mix1ph->Pz(),mix1ph->E());
587 Double_t pi0gammapt=(vpi0+vmixph).Pt();
588 Double_t pi0gammamass=(vpi0+vmixph).M();
589 Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
590 Double_t gammaOverOmegaPtRatio= vmixph.Pt()/pi0gammapt;
591
592 //pi0, gamma pt cut
593 if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
594 gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
595
596 for(Int_t ipid=0;ipid<fNpid;ipid++){
597 for(Int_t idist=0;idist<fNBadChDistBin;idist++){
598 Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
599 if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
600 photon2->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
601 mix1ph->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
602 photon1->DistToBad()>=idist &&
603 photon2->DistToBad()>=idist &&
604 mix1ph->DistToBad()>=idist ){
605 if(GetDebug() > 2) printf("MixA: index %d pt %2.3f mass %2.3f \n",index, pi0gammapt, pi0gammamass);
606 //fill the histograms
607 fMixAOmega0[index]->Fill(pi0gammapt,pi0gammamass);
608 if(pi0asy<0.7)fMixAOmega1[index]->Fill(pi0gammapt,pi0gammamass);
609 if(pi0asy<0.8)fMixAOmega2[index]->Fill(pi0gammapt,pi0gammamass);
610 //printf("mix A %d %2.2f \n", index, pi0gammamass);
611
612 }
613 }
614 }
615 }
616 }
617 }
618 }
619
620 //
621 // --B (r1_event1+r2_event2)+r3_event2
622 //
623 if(GetDebug() >0)printf("MixB: (r1_event1+r2_event2)+r3_event2 \n");
624 for(Int_t i=0;i<nphotons;i++){
625 AliAODPWG4Particle *ph1 = (AliAODPWG4Particle*) (aodGamma->At(i));
626 TLorentzVector vph1(ph1->Px(),ph1->Py(),ph1->Pz(),ph1->E());
627 //interrupt here...................
628 //especially for EMCAL
629 //we suppose the high pt clusters are overlapped pi0
630
631 for(Int_t j=i+1;j<nphotons;j++){
632 AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (aodGamma->At(j));
633 TLorentzVector fakePi0, fakeOmega, vph;
634
635 if(ph1->E() > fEOverlapCluster && ph1->E() > ph2->E()) {
636 fakePi0.SetPxPyPzE(ph1->Px(),ph1->Py(),ph1->Pz(),TMath::Sqrt(ph1->Px()*ph1->Px()+ph1->Py()*ph1->Py()+ph1->Pz()*ph1->Pz()+0.135*0.135));
637 vph.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
638 }
639 else if(ph2->E() > fEOverlapCluster && ph2->E() > ph1->E()) {
640 fakePi0.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),TMath::Sqrt(ph2->Px()*ph2->Px()+ph2->Py()*ph2->Py()+ph2->Pz()*ph2->Pz()+0.135*0.135));
641 vph.SetPxPyPzE(ph1->Px(), ph1->Py(),ph1->Pz(),ph1->E());
642 }
643 else continue;
644
645 fakeOmega=fakePi0+vph;
646 for(Int_t ii=0;ii<fNCentBin;ii++){
647 fhFakeOmega[icentbin]->Fill(fakeOmega.Pt(), fakeOmega.M());
648 }
649 }//j
650
651 //continue ................
652 Int_t nMixed = fEventsList[curEventBin]->GetSize();
653 for(Int_t ie=0;ie<nMixed;ie++){
654 TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(ie));
655 for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
656 AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (ev2->At(mix1));
657 TLorentzVector vph2(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
658 Double_t pi0asy = TMath::Abs(vph1.E()-vph2.E())/(vph1.E()+vph2.E());
659 Double_t pi0mass=(vph1+vph2).M();
660
661 if(TMath::Abs(pi0mass-fPi0Mass)<fPi0MassWindow){//for pi0 selection
662 for(Int_t mix2=(mix1+1);mix2<ev2->GetEntries();mix2++){
663 AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev2->At(mix2));
664 TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
665
666 Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
667 Double_t pi0gammamass=(vph1+vph2+vph3).M();
668 Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
669 Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
670 //pi0, gamma pt cut
671 if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
672 gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
673
674 for(Int_t ipid=0;ipid<fNpid;ipid++){
675 for(Int_t idist=0;idist<fNBadChDistBin;idist++){
676 Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
677 if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
678 ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
679 ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
680 ph1->DistToBad()>=idist &&
681 ph2->DistToBad()>=idist &&
682 ph3->DistToBad()>=idist ){
683 if(GetDebug() > 2) printf("MixB: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
684 //fill histograms
685 fMixBOmega0[index]->Fill(pi0gammapt,pi0gammamass);
686 if(pi0asy<0.7) fMixBOmega1[index]->Fill(pi0gammapt,pi0gammamass);
687 if(pi0asy<0.8) fMixBOmega2[index]->Fill(pi0gammapt,pi0gammamass);
688 //printf("mix B %d %2.2f \n", index, pi0gammamass);
689 }
690 }
691 }
692 }
693
694 //
695 // --C (r1_event1+r2_event2)+r3_event3
696 //
697 if(GetDebug() >0)printf("MixC: (r1_event1+r2_event2)+r3_event3\n");
698 for(Int_t je=(ie+1);je<nMixed;je++){
699 TClonesArray* ev3= (TClonesArray*) (fEventsList[curEventBin]->At(je));
700 for(Int_t mix3=0;mix3<ev3->GetEntries();mix3++){
701 AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev3->At(mix3));
702 TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
703
704 Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
705 Double_t pi0gammamass=(vph1+vph2+vph3).M();
706 Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
707 Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
708 //pi0, gamma pt cut
709 if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
710 gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
711
712 for(Int_t ipid=0;ipid<fNpid;ipid++){
713 for(Int_t idist=0;idist<fNBadChDistBin;idist++){
714 Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
715 if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
716 ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
717 ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
718 ph1->DistToBad()>=idist &&
719 ph2->DistToBad()>=idist &&
720 ph3->DistToBad()>=idist ){
721 if(GetDebug() > 2) printf("MixC: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
722 //fill histograms
723 fMixCOmega0[index]->Fill(pi0gammapt,pi0gammamass);
724 if(pi0asy<0.7) fMixCOmega1[index]->Fill(pi0gammapt,pi0gammamass);
725 if(pi0asy<0.8) fMixCOmega2[index]->Fill(pi0gammapt,pi0gammamass);
726 //printf("mix C %d %2.2f \n", index, pi0gammamass);
727 }
728 }
729 }
730 }
731 }
732 } //for pi0 selecton
733 }
734 }
735 }
736
737
738 //event buffer
739 TClonesArray *currentEvent = new TClonesArray(*aodGamma);
740 if(currentEvent->GetEntriesFast()>0){
741 fEventsList[curEventBin]->AddFirst(currentEvent) ;
742 currentEvent=0 ;
743 if(fEventsList[curEventBin]->GetSize()>=fNmaxMixEv) {
744 TClonesArray * tmp = (TClonesArray*) (fEventsList[curEventBin]->Last()) ;
745 fEventsList[curEventBin]->RemoveLast() ;
746 delete tmp ;
747 }
748 }
749 else{
750 delete currentEvent ;
751 currentEvent=0 ;
752 }
753
754}
755
756//______________________________________________________________________________
757void AliAnaOmegaToPi0Gamma::ReadHistograms(TList * outputList)
758{
759 //read the histograms
760 //for the finalization of the terminate analysis
761
762 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"RealToPi0Gamma_Vz0C0Rp0Pid0Dist0"));
763
764 Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
765
766 if(!fRealOmega0) fRealOmega0 =new TH2F*[ndim];
767 if(!fMixAOmega0) fMixAOmega0 =new TH2F*[ndim];
768 if(!fMixBOmega0) fMixBOmega0 =new TH2F*[ndim];
769 if(!fMixCOmega0) fMixCOmega0 =new TH2F*[ndim];
770
771 if(!fRealOmega1) fRealOmega1 =new TH2F*[ndim];
772 if(!fMixAOmega1) fMixAOmega1 =new TH2F*[ndim];
773 if(!fMixBOmega1) fMixBOmega1 =new TH2F*[ndim];
774 if(!fMixCOmega1) fMixCOmega1 =new TH2F*[ndim];
775
776 if(!fRealOmega2) fRealOmega2 =new TH2F*[ndim];
777 if(!fMixAOmega2) fMixAOmega2 =new TH2F*[ndim];
778 if(!fMixBOmega2) fMixBOmega2 =new TH2F*[ndim];
779 if(!fMixCOmega2) fMixCOmega2 =new TH2F*[ndim];
780
781 for(Int_t i=0;i<fNVtxZBin;i++){
782 for(Int_t j=0;j<fNCentBin;j++){
783 for(Int_t k=0;k<fNRpBin;k++){ //at event level
784 Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
785 for(Int_t ipid=0;ipid<fNpid;ipid++){
786 for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle
787 Int_t ind=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
788 fRealOmega0[ind]= (TH2F*) outputList->At(index++);
789 fMixAOmega0[ind]= (TH2F*) outputList->At(index++);
790 fMixBOmega0[ind]= (TH2F*) outputList->At(index++);
791 fMixCOmega0[ind]= (TH2F*) outputList->At(index++);
792
793 fRealOmega1[ind]= (TH2F*) outputList->At(index++);
794 fMixAOmega1[ind]= (TH2F*) outputList->At(index++);
795 fMixBOmega1[ind]= (TH2F*) outputList->At(index++);
796 fMixCOmega1[ind]= (TH2F*) outputList->At(index++);
797
798 fRealOmega2[ind]= (TH2F*) outputList->At(index++);
799 fMixAOmega2[ind]= (TH2F*) outputList->At(index++);
800 fMixBOmega2[ind]= (TH2F*) outputList->At(index++);
801 fMixCOmega2[ind]= (TH2F*) outputList->At(index++);
802
803
804 }
805 }
806 }
807 }
808 }
809
810 if(IsDataMC()){
811 fhOmegaPriPt = (TH1F*) outputList->At(index++);
812 }
813
814}
815
816//______________________________________________________________________________
817void AliAnaOmegaToPi0Gamma::Terminate(TList * outputList)
818{
819// //Do some calculations and plots from the final histograms.
820 if(GetDebug() >= 0) printf("AliAnaOmegaToPi0Gamma::Terminate() \n");
821 ReadHistograms(outputList);
822 const Int_t buffersize = 255;
823 char cvs1[buffersize];
824 snprintf(cvs1, buffersize, "Neutral_%s_IVM",fInputAODGammaName.Data());
825
826 TCanvas * cvsIVM = new TCanvas(cvs1, cvs1, 400, 10, 600, 700) ;
827 cvsIVM->Divide(2, 2);
828
829 cvsIVM->cd(1);
830 char dec[buffersize];
831 snprintf(dec,buffersize,"h2Real_%s",fInputAODGammaName.Data());
832 TH2F * h2Real= (TH2F*)fRealOmega0[0]->Clone(dec);
833 h2Real->GetXaxis()->SetRangeUser(4,6);
834 TH1F * hRealOmega = (TH1F*) h2Real->ProjectionY();
835 hRealOmega->SetTitle("RealPi0Gamma 4<pt<6");
836 hRealOmega->SetLineColor(2);
837 hRealOmega->Draw();
838
839 cvsIVM->cd(2);
840 snprintf(dec,buffersize,"hMixA_%s",fInputAODGammaName.Data());
841 TH2F *h2MixA= (TH2F*)fMixAOmega0[0]->Clone(dec);
842 h2MixA->GetXaxis()->SetRangeUser(4,6);
843 TH1F * hMixAOmega = (TH1F*) h2MixA->ProjectionY();
844 hMixAOmega->SetTitle("MixA 4<pt<6");
845 hMixAOmega->SetLineColor(2);
846 hMixAOmega->Draw();
847
848 cvsIVM->cd(3);
849 snprintf(dec,buffersize,"hMixB_%s",fInputAODGammaName.Data());
850 TH2F * h2MixB= (TH2F*)fMixBOmega0[0]->Clone(dec);
851 h2MixB->GetXaxis()->SetRangeUser(4,6);
852 TH1F * hMixBOmega = (TH1F*) h2MixB->ProjectionY();
853 hMixBOmega->SetTitle("MixB 4<pt<6");
854 hMixBOmega->SetLineColor(2);
855 hMixBOmega->Draw();
856
857 cvsIVM->cd(4);
858 snprintf(dec,buffersize,"hMixC_%s",fInputAODGammaName.Data());
859 TH2F *h2MixC= (TH2F*)fMixCOmega0[0]->Clone(dec);
860 h2MixC->GetXaxis()->SetRangeUser(4,6);
861 TH1F * hMixCOmega = (TH1F*) h2MixC->ProjectionY();
862 hMixCOmega->SetTitle("MixC 4<pt<6");
863 hMixCOmega->SetLineColor(2);
864 hMixCOmega->Draw();
865
866 char eps[buffersize];
867 snprintf(eps,buffersize,"CVS_%s_IVM.eps",fInputAODGammaName.Data());
868 cvsIVM->Print(eps);
869 cvsIVM->Modified();
870
871}