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