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