AliPhysics  e34b7ac (e34b7ac)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaOmegaToPi0Gamma.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system
17 class TROOT;
18 #include "TH2F.h"
19 #include "TLorentzVector.h"
20 #include "TParticle.h"
21 #include "TCanvas.h"
22 #include "TFile.h"
23 
24 //---- AliRoot system ----
25 #include "AliAnaOmegaToPi0Gamma.h"
26 #include "AliCaloTrackReader.h"
27 #include "AliCaloPID.h"
28 #include "AliStack.h"
29 #include "AliVEvent.h"
30 #include "AliAODEvent.h"
31 #include "AliAODMCParticle.h"
32 
36 
37 //______________________________________________________________________________
39 //______________________________________________________________________________
41 fInputAODPi0(0), fInputAODGammaName(""),
42 fEventsList(0x0),fNVtxZBin(0), fNCentBin(0), fNRpBin(0), fNBadChDistBin(0), fNpid(0),
43 fVtxZCut(0), fCent(0), fRp(0),
44 fPi0Mass(0),fPi0MassWindow(0),fPi0OverOmegaPtCut(0),
45 fGammaOverOmegaPtCut(0), fEOverlapCluster(0),
46 fhEtalon(0),
47 fRealOmega0(0), fMixAOmega0(0),
48 fMixBOmega0(0), fMixCOmega0(0),
49 fRealOmega1(0), fMixAOmega1(0),
50 fMixBOmega1(0), fMixCOmega1(0),
51 fRealOmega2(0), fMixAOmega2(0),
52 fMixBOmega2(0), fMixCOmega2(0),
53 fhFakeOmega(0),
54 fhOmegaPriPt(0)
55 {
57 }
58 
59 //_____________________________________________
61 //_____________________________________________
63 {
64 // Done by the maker
65 // if(fInputAODPi0){
66 // fInputAODPi0->Clear();
67 // delete fInputAODPi0;
68 // }
69 
70  if(fEventsList){
71  for(Int_t i=0;i<fNVtxZBin;i++)
72  {
73  for(Int_t j=0;j<fNCentBin;j++)
74  {
75  for(Int_t k=0;k<fNRpBin;k++)
76  {
77  fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->Clear();
78  delete fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k];
79  }
80  }
81  }
82  }
83 
84  delete [] fEventsList;
85  fEventsList=0;
86 
87  delete [] fVtxZCut;
88  delete [] fCent;
89  delete [] fRp;
90 }
91 
92 //__________________________________________
95 //__________________________________________
97 {
98  fInputAODGammaName="PhotonsDetector";
99  fNVtxZBin=1;
100  fNCentBin=1;
101  fNRpBin=1;
102  fNBadChDistBin=3;
103  fNpid=1;
104 
105  fPi0Mass=0.1348;
106  fPi0MassWindow=0.015;
107  fPi0OverOmegaPtCut=0.8;
109 
111 }
112 
113 //_____________________________________________________
115 //_____________________________________________________
117 {
118  fVtxZCut = new Double_t [fNVtxZBin];
119  for(Int_t i=0;i<fNVtxZBin;i++) fVtxZCut[i]=10*(i+1);
120 
121  fCent=new Double_t[fNCentBin];
122  for(int i = 0;i<fNCentBin;i++)fCent[i]=0;
123 
124  fRp=new Double_t[fNRpBin];
125  for(int i = 0;i<fNRpBin;i++)fRp[i]=0;
126  //
130 
131  Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
134 
135  fhEtalon = new TH2F("hEtalon","Histo with binning parameters", nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
136  fhEtalon->SetXTitle("P_{T} (GeV)") ;
137  fhEtalon->SetYTitle("m_{inv} (GeV)") ;
138 
139  // store them in fOutputContainer
140  fEventsList = new TList*[fNVtxZBin*fNCentBin*fNRpBin];
141  for(Int_t i=0;i<fNVtxZBin;i++){
142  for(Int_t j=0;j<fNCentBin;j++){
143  for(Int_t k=0;k<fNRpBin;k++){
144  fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k] = new TList();
145  fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->SetOwner(kFALSE);
146  }
147  }
148  }
149 
150  TList * outputContainer = new TList() ;
151  outputContainer->SetName(GetName());
152  const Int_t buffersize = 255;
153  char key[buffersize] ;
154  char title[buffersize] ;
155  const char * detectorName= fInputAODGammaName.Data();
156  Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
157 
158  fRealOmega0 =new TH2F*[ndim];
159  fMixAOmega0 =new TH2F*[ndim];
160  fMixBOmega0 =new TH2F*[ndim];
161  fMixCOmega0 =new TH2F*[ndim];
162 
163  fRealOmega1 =new TH2F*[ndim];
164  fMixAOmega1 =new TH2F*[ndim];
165  fMixBOmega1 =new TH2F*[ndim];
166  fMixCOmega1 =new TH2F*[ndim];
167 
168  fRealOmega2 =new TH2F*[ndim];
169  fMixAOmega2 =new TH2F*[ndim];
170  fMixBOmega2 =new TH2F*[ndim];
171  fMixCOmega2 =new TH2F*[ndim];
172 
173  fhFakeOmega = new TH2F*[fNCentBin];
174 
175  for(Int_t i=0;i<fNVtxZBin;i++){
176  for(Int_t j=0;j<fNCentBin;j++){
177  for(Int_t k=0;k<fNRpBin;k++){ //at event level
178  Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
179  for(Int_t ipid=0;ipid<fNpid;ipid++){
180  for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle level
181 
182  Int_t index=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
183 
184  snprintf(key,buffersize,"RealToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
185  snprintf(title,buffersize, "%s Real Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
186  fRealOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
187  fRealOmega0[index]->SetName(key) ;
188  fRealOmega0[index]->SetTitle(title);
189  outputContainer->Add(fRealOmega0[index]);
190 
191  snprintf(key,buffersize,"MixAToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
192  snprintf(title,buffersize, "%s MixA Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
193  fMixAOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
194  fMixAOmega0[index]->SetName(key) ;
195  fMixAOmega0[index]->SetTitle(title);
196  outputContainer->Add(fMixAOmega0[index]);
197 
198  snprintf(key,buffersize,"MixBToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
199  snprintf(title,buffersize, "%s MixB Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
200  fMixBOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
201  fMixBOmega0[index]->SetName(key) ;
202  fMixBOmega0[index]->SetTitle(title);
203  outputContainer->Add(fMixBOmega0[index]);
204 
205  snprintf(key,buffersize,"MixCToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
206  snprintf(title,buffersize, "%s MixC Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
207  fMixCOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
208  fMixCOmega0[index]->SetName(key) ;
209  fMixCOmega0[index]->SetTitle(title);
210  outputContainer->Add(fMixCOmega0[index]);
211 
212  snprintf(key,buffersize,"RealToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
213  snprintf(title,buffersize, "%s Real Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
214  fRealOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
215  fRealOmega1[index]->SetName(key) ;
216  fRealOmega1[index]->SetTitle(title);
217  outputContainer->Add(fRealOmega1[index]);
218 
219  snprintf(key,buffersize,"MixAToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
220  snprintf(title,buffersize, "%s MixA Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
221  fMixAOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
222  fMixAOmega1[index]->SetName(key) ;
223  fMixAOmega1[index]->SetTitle(title);
224  outputContainer->Add(fMixAOmega1[index]);
225 
226  snprintf(key,buffersize,"MixBToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
227  snprintf(title,buffersize, "%s MixB Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
228  fMixBOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
229  fMixBOmega1[index]->SetName(key) ;
230  fMixBOmega1[index]->SetTitle(title);
231  outputContainer->Add(fMixBOmega1[index]);
232 
233  snprintf(key,buffersize,"MixCToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
234  snprintf(title,buffersize, "%s MixC Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
235  fMixCOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
236  fMixCOmega1[index]->SetName(key) ;
237  fMixCOmega1[index]->SetTitle(title);
238  outputContainer->Add(fMixCOmega1[index]);
239 
240  snprintf(key,buffersize,"RealToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
241  snprintf(title,buffersize, "%s Real Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
242  fRealOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
243  fRealOmega2[index]->SetName(key) ;
244  fRealOmega2[index]->SetTitle(title);
245  outputContainer->Add(fRealOmega2[index]);
246 
247  snprintf(key,buffersize,"MixAToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
248  snprintf(title,buffersize, "%s MixA Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
249  fMixAOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
250  fMixAOmega2[index]->SetName(key) ;
251  fMixAOmega2[index]->SetTitle(title);
252  outputContainer->Add(fMixAOmega2[index]);
253 
254  snprintf(key,buffersize,"MixBToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
255  snprintf(title,buffersize, "%s MixB Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
256  fMixBOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
257  fMixBOmega2[index]->SetName(key) ;
258  fMixBOmega2[index]->SetTitle(title);
259  outputContainer->Add(fMixBOmega2[index]);
260 
261  snprintf(key,buffersize,"MixCToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
262  snprintf(title,buffersize, "%s MixC Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detectorName,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
263  fMixCOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
264  fMixCOmega2[index]->SetName(key) ;
265  fMixCOmega2[index]->SetTitle(title);
266  outputContainer->Add(fMixCOmega2[index]);
267  }
268  }
269  }
270  }
271  }
272 
273  for(Int_t i=0;i<fNCentBin;i++){
274  snprintf(key,buffersize, "fhFakeOmega%d",i);
275  snprintf(title,buffersize,"FakePi0(high pt cluster)+Gamma with Centrality%d",i);
276  fhFakeOmega[i] = (TH2F*)fhEtalon->Clone(key);
277  fhFakeOmega[i]->SetTitle(title);
278  outputContainer->Add(fhFakeOmega[i]);
279  }
280  if(IsDataMC()){
281  snprintf(key,buffersize, "%sOmegaPri",detectorName);
282  snprintf(title,buffersize,"primary #omega in %s",detectorName);
283  fhOmegaPriPt=new TH1F(key, title,nptbins,ptmin,ptmax);
284  fhOmegaPriPt->GetXaxis()->SetTitle("P_{T}");
285  fhOmegaPriPt->GetYaxis()->SetTitle("dN/P_{T}");
286  outputContainer->Add(fhOmegaPriPt);
287  }
288 
289  delete fhEtalon;
290  return outputContainer;
291 }
292 
293 //_______________________________________________________________
295 //_______________________________________________________________
296 void AliAnaOmegaToPi0Gamma::Print(const Option_t * /*opt*/) const
297 {
298  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
300  printf("Omega->pi0+gamma->3gamma\n");
301  printf("Cuts at event level: \n");
302  printf("Bins of vertex Z: %d \n", fNVtxZBin);
303  printf("Bins of centrality: %d \n",fNCentBin);
304  printf("Bins of Reaction plane: %d\n",fNRpBin);
305  printf("Cuts at AOD particle level:\n");
306  printf("Number of PID: %d \n", fNpid);
307  printf("Number of DistToBadChannel cuts: %d\n", fNBadChDistBin);
308 }
309 
310 //______________________________________________________
312 //______________________________________________________
314 {
315  // Fill the MC AOD if needed first.
316  //-----------
317  //need to be further implemented
318  AliStack * stack = 0x0;
319  // TParticle * primary = 0x0;
320  TClonesArray * mcparticles0 = 0x0;
321  //TClonesArray * mcparticles1 = 0x0;
322  AliAODMCParticle * aodprimary = 0x0;
323  Int_t pdg=0;
324  Double_t pt=0;
325  Double_t eta=0;
326 
327  if(IsDataMC())
328  {
329  if(GetReader()->ReadStack())
330  {
331  stack = GetMCStack() ;
332  if(!stack){
333  printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - There is no stack!\n");
334  }
335  else{
336  for(Int_t i=0 ; i<stack->GetNtrack(); i++){
337  TParticle * prim = stack->Particle(i) ;
338  pdg = prim->GetPdgCode() ;
339  eta=prim->Eta();
340  pt=prim->Pt();
341  if(TMath::Abs(eta)<0.5) {
342  if(pdg==223) fhOmegaPriPt->Fill(pt, GetEventWeight());
343  }
344  }
345  }
346  }
347  else if(GetReader()->ReadAODMCParticles()){
348  //Get the list of MC particles
349  mcparticles0 = GetReader()->GetAODMCParticles();
350  if(!mcparticles0 )
351  {
352  if(GetDebug() > 0) printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
353  }
354  else
355  {
356  for(Int_t i=0;i<mcparticles0->GetEntries();i++)
357  {
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  {
364  if(pdg==223) fhOmegaPriPt->Fill(pt, GetEventWeight());
365  }
366 
367  }
368  }//mcparticles0 exists
369  }//AOD MC Particles
370  }// is data and MC
371 
372 
373  //process event from AOD brach
374  //extract pi0, eta and omega analysis
375  Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
376  if(IsBadRun(iRun)) return ;
377 
378  //vertex z
379  Double_t vert[]={0,0,0} ;
380  GetVertex(vert);
381  Int_t curEventBin =0;
382 
383  Int_t ivtxzbin=(Int_t)TMath::Abs(vert[2])/10;
384  if(ivtxzbin>=fNVtxZBin)return;
385 
386  //centrality
387  Int_t currentCentrality = GetEventCentrality();
388  if(currentCentrality == -1) return;
389  Int_t optCent = GetReader()->GetCentralityOpt();
390  Int_t icentbin=currentCentrality/(optCent/fNCentBin) ; //GetEventCentrality();
391 
392  printf("-------------- %d %d %d ",currentCentrality, optCent, icentbin);
393  //reaction plane
394  Int_t irpbin=0;
395 
396  if(ivtxzbin==-1) return;
397  curEventBin = ivtxzbin*fNCentBin*fNRpBin + icentbin*fNRpBin + irpbin;
398  TClonesArray *aodGamma = (TClonesArray*) GetAODBranch(fInputAODGammaName); //photon array
399  // TClonesArray *aodGamma = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaName); //photon array
400  Int_t nphotons =0;
401  if(aodGamma) nphotons= aodGamma->GetEntries();
402  else return;
403 
404  fInputAODPi0 = (TClonesArray*)GetInputAODBranch(); //pi0 array
405  Int_t npi0s = 0;
406  if(fInputAODPi0) npi0s= fInputAODPi0 ->GetEntries();
407  else return;
408 
409  if(nphotons<3 || npi0s<1)return; //for pi0, eta and omega->pi0+gamma->3gamma reconstruction
410 
411  //reconstruction of omega(782)->pi0+gamma->3gamma
412  //loop for pi0 and photon
413  if(GetDebug() > 0) printf("omega->pi0+gamma->3gamma invariant mass analysis ! This event have %d photons and %d pi0 \n", nphotons, npi0s);
414  for(Int_t i=0;i<npi0s;i++){
415  AliAODPWG4Particle * pi0 = (AliAODPWG4Particle*) (fInputAODPi0->At(i)) ; //pi0
416  TLorentzVector vpi0(pi0->Px(),pi0->Py(),pi0->Pz(),pi0->E());
417  Int_t lab1=pi0->GetCaloLabel(0); // photon1 from pi0 decay
418  Int_t lab2=pi0->GetCaloLabel(1); // photon2 from pi0 decay
419  //for omega->pi0+gamma, it needs at least three photons per event
420  //Get the two decay photons from pi0
421  AliAODPWG4Particle * photon1 =0;
422  AliAODPWG4Particle * photon2 =0;
423  for(Int_t d1=0;d1<nphotons;d1++){
424  for(Int_t d2=0;d2<nphotons;d2++){
425  AliAODPWG4Particle * dp1 = (AliAODPWG4Particle*) (aodGamma->At(d1));
426  AliAODPWG4Particle * dp2 = (AliAODPWG4Particle*) (aodGamma->At(d2));
427  Int_t dlab1=dp1->GetCaloLabel(0);
428  Int_t dlab2=dp2->GetCaloLabel(0);
429  if(dlab1==lab1 && dlab2==lab2){
430  photon1=dp1;
431  photon2=dp2;
432  }
433  else continue;
434  }
435  }
436  //caculate the asy and dist of the two photon from pi0 decay
437  TLorentzVector dph1(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());
438  TLorentzVector dph2(photon2->Px(),photon2->Py(),photon2->Pz(),photon2->E());
439 
440  Double_t pi0asy= TMath::Abs(dph1.E()-dph2.E())/(dph1.E()+dph2.E());
441  // Double_t phi1=dph1.Phi();
442  // Double_t phi2=dph2.Phi();
443  // Double_t eta1=dph1.Eta();
444  // Double_t eta2=dph2.Eta();
445  // Double_t pi0dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));
446 
447  if(pi0->GetIdentifiedParticleType()==111 && nphotons>2 && npi0s
448  && TMath::Abs(vpi0.M()-fPi0Mass)<fPi0MassWindow) { //pi0 candidates
449 
450  //avoid the double counting
451  Int_t * dc1= new Int_t[nphotons];
452  Int_t * dc2= new Int_t[nphotons];
453  Int_t index1=0;
454  Int_t index2=0;
455  for(Int_t k=0;k<i;k++){
456  AliAODPWG4Particle * p3=(AliAODPWG4Particle*)(fInputAODPi0->At(k));
457  Int_t lab4=p3->GetCaloLabel(0);
458  Int_t lab5=p3->GetCaloLabel(1);
459  if(lab1==lab4){ dc1[index1]=lab5; index1++; }
460  if(lab2==lab5){ dc2[index2]=lab4; index2++; }
461  }
462 
463 
464  //loop the pi0 with third gamma
465  for(Int_t j=0;j<nphotons;j++){
466  AliAODPWG4Particle *photon3 = (AliAODPWG4Particle*) (aodGamma->At(j));
467  TLorentzVector dph3(photon3->Px(),photon3->Py(),photon3->Pz(),photon3->E());
468  Int_t lab3=photon3->GetCaloLabel(0);
469  Double_t pi0gammapt=(vpi0+dph3).Pt();
470  Double_t pi0gammamass=(vpi0+dph3).M();
471  Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
472  Double_t gammaOverOmegaPtRatio= dph3.Pt()/pi0gammapt;
473 
474  //pi0, gamma pt cut
475  if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
476  gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
477 
478  for(Int_t l=0;l<index1;l++) if(lab3==dc1[l]) lab3=-1;
479  for(Int_t l=0;l<index2;l++) if(lab3==dc2[l]) lab3=-1;
480 
481  if(lab3>0 && lab3!=lab1 && lab3!=lab2){
482  for(Int_t ipid=0;ipid<fNpid;ipid++){
483  for(Int_t idist=0;idist<fNBadChDistBin;idist++){
484  Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
485  if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
486  photon2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
487  photon3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
488  photon1->DistToBad()>=idist &&
489  photon2->DistToBad()>=idist &&
490  photon3->DistToBad()>=idist ){
491  //fill the histograms
492  if(GetDebug() > 2) printf("Real: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
493  fRealOmega0[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
494  if(pi0asy<0.7) fRealOmega1[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
495  if(pi0asy<0.8) fRealOmega2[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
496  }
497  }
498  }
499  }
500  }
501  delete []dc1;
502  delete []dc2;
503  if(GetDebug() > 0) printf("MixA: (r1_event1+r2_event1)+r3_event2 \n");
504  //-------------------------
505  //background analysis
506  //three background
507  // --A (r1_event1+r2_event1)+r3_event2
508  Int_t nMixed = fEventsList[curEventBin]->GetSize();
509  for(Int_t im=0;im<nMixed;im++){
510  TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(im));
511  for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
512  AliAODPWG4Particle *mix1ph = (AliAODPWG4Particle*) (ev2->At(mix1));
513  TLorentzVector vmixph(mix1ph->Px(),mix1ph->Py(),mix1ph->Pz(),mix1ph->E());
514  Double_t pi0gammapt=(vpi0+vmixph).Pt();
515  Double_t pi0gammamass=(vpi0+vmixph).M();
516  Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
517  Double_t gammaOverOmegaPtRatio= vmixph.Pt()/pi0gammapt;
518 
519  //pi0, gamma pt cut
520  if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
521  gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
522 
523  for(Int_t ipid=0;ipid<fNpid;ipid++){
524  for(Int_t idist=0;idist<fNBadChDistBin;idist++){
525  Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
526  if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
527  photon2->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
528  mix1ph->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
529  photon1->DistToBad()>=idist &&
530  photon2->DistToBad()>=idist &&
531  mix1ph->DistToBad()>=idist ){
532  if(GetDebug() > 2) printf("MixA: index %d pt %2.3f mass %2.3f \n",index, pi0gammapt, pi0gammamass);
533  //fill the histograms
534  fMixAOmega0[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
535  if(pi0asy<0.7)fMixAOmega1[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
536  if(pi0asy<0.8)fMixAOmega2[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
537  //printf("mix A %d %2.2f \n", index, pi0gammamass);
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(), GetEventWeight());
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, GetEventWeight());
612  if(pi0asy<0.7) fMixBOmega1[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
613  if(pi0asy<0.8) fMixBOmega2[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
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, GetEventWeight());
650  if(pi0asy<0.7) fMixCOmega1[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
651  if(pi0asy<0.8) fMixCOmega2[index]->Fill(pi0gammapt, pi0gammamass, GetEventWeight());
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 ;
669  if(fEventsList[curEventBin]->GetSize()>=GetNMaxEvMix()) {
670  TClonesArray * tmp = (TClonesArray*) (fEventsList[curEventBin]->Last()) ;
671  fEventsList[curEventBin]->RemoveLast() ;
672  delete tmp ;
673  }
674  }
675  else
676  {
677  delete currentEvent ;
678  currentEvent=0 ;
679  }
680 }
681 
682 //____________________________________________________________
685 //____________________________________________________________
687 {
688  Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"RealToPi0Gamma_Vz0C0Rp0Pid0Dist0"));
689 
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  {
709  for(Int_t j=0;j<fNCentBin;j++)
710  {
711  for(Int_t k=0;k<fNRpBin;k++)
712  { //at event level
713  Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
714 
715  for(Int_t ipid=0;ipid<fNpid;ipid++)
716  {
717  for(Int_t idist=0;idist<fNBadChDistBin;idist++)
718  { //at particle
719  Int_t ind=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
720 
721  fRealOmega0[ind]= (TH2F*) outputList->At(index++);
722  fMixAOmega0[ind]= (TH2F*) outputList->At(index++);
723  fMixBOmega0[ind]= (TH2F*) outputList->At(index++);
724  fMixCOmega0[ind]= (TH2F*) outputList->At(index++);
725 
726  fRealOmega1[ind]= (TH2F*) outputList->At(index++);
727  fMixAOmega1[ind]= (TH2F*) outputList->At(index++);
728  fMixBOmega1[ind]= (TH2F*) outputList->At(index++);
729  fMixCOmega1[ind]= (TH2F*) outputList->At(index++);
730 
731  fRealOmega2[ind]= (TH2F*) outputList->At(index++);
732  fMixAOmega2[ind]= (TH2F*) outputList->At(index++);
733  fMixBOmega2[ind]= (TH2F*) outputList->At(index++);
734  fMixCOmega2[ind]= (TH2F*) outputList->At(index++);
735  }
736  }
737  }
738  }
739  }
740 
741  if(IsDataMC())
742  {
743  fhOmegaPriPt = (TH1F*) outputList->At(index++);
744  }
745 }
746 
747 //_______________________________________________________
749 //_______________________________________________________
751 {
752  if(GetDebug() >= 0) printf("AliAnaOmegaToPi0Gamma::Terminate() \n");
753  ReadHistograms(outputList);
754  const Int_t buffersize = 255;
755  char cvs1[buffersize];
756  snprintf(cvs1, buffersize, "Neutral_%s_IVM",fInputAODGammaName.Data());
757 
758  TCanvas * cvsIVM = new TCanvas(cvs1, cvs1, 400, 10, 600, 700) ;
759  cvsIVM->Divide(2, 2);
760 
761  cvsIVM->cd(1);
762  char dec[buffersize];
763  snprintf(dec,buffersize,"h2Real_%s",fInputAODGammaName.Data());
764  TH2F * h2Real= (TH2F*)fRealOmega0[0]->Clone(dec);
765  h2Real->GetXaxis()->SetRangeUser(4,6);
766  TH1F * hRealOmega = (TH1F*) h2Real->ProjectionY();
767  hRealOmega->SetTitle("RealPi0Gamma 4<pt<6");
768  hRealOmega->SetLineColor(2);
769  hRealOmega->Draw();
770 
771  cvsIVM->cd(2);
772  snprintf(dec,buffersize,"hMixA_%s",fInputAODGammaName.Data());
773  TH2F *h2MixA= (TH2F*)fMixAOmega0[0]->Clone(dec);
774  h2MixA->GetXaxis()->SetRangeUser(4,6);
775  TH1F * hMixAOmega = (TH1F*) h2MixA->ProjectionY();
776  hMixAOmega->SetTitle("MixA 4<pt<6");
777  hMixAOmega->SetLineColor(2);
778  hMixAOmega->Draw();
779 
780  cvsIVM->cd(3);
781  snprintf(dec,buffersize,"hMixB_%s",fInputAODGammaName.Data());
782  TH2F * h2MixB= (TH2F*)fMixBOmega0[0]->Clone(dec);
783  h2MixB->GetXaxis()->SetRangeUser(4,6);
784  TH1F * hMixBOmega = (TH1F*) h2MixB->ProjectionY();
785  hMixBOmega->SetTitle("MixB 4<pt<6");
786  hMixBOmega->SetLineColor(2);
787  hMixBOmega->Draw();
788 
789  cvsIVM->cd(4);
790  snprintf(dec,buffersize,"hMixC_%s",fInputAODGammaName.Data());
791  TH2F *h2MixC= (TH2F*)fMixCOmega0[0]->Clone(dec);
792  h2MixC->GetXaxis()->SetRangeUser(4,6);
793  TH1F * hMixCOmega = (TH1F*) h2MixC->ProjectionY();
794  hMixCOmega->SetTitle("MixC 4<pt<6");
795  hMixCOmega->SetLineColor(2);
796  hMixCOmega->Draw();
797 
798  char eps[buffersize];
799  snprintf(eps,buffersize,"CVS_%s_IVM.eps",fInputAODGammaName.Data());
800  cvsIVM->Print(eps);
801  cvsIVM->Modified();
802 }
Float_t GetHistoPtMax() const
TH2F ** fRealOmega1
! Real omega IVM(asy, pt, m), with Asy_pi0<0.7
Int_t pdg
TH2F ** fMixCOmega0
! MixC omega IVM(asy, pt, m)
Float_t GetHistoPtMin() const
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
virtual void GetVertex(Double_t vertex[3]) const
virtual AliVEvent * GetInputEvent() const
TList * GetCreateOutputObjects()
Create histograms to be saved in output file.
Int_t fNRpBin
Number of reaction plane cut.
TH2F ** fhFakeOmega
! High pt clusters assumed as pi0 + another gamma
void Terminate(TList *outList)
Do some calculations and plots from the final histograms.
Double_t fGammaOverOmegaPtCut
gamma pt over omega pt cut
Int_t GetHistoMassBins() const
Float_t GetHistoMassMin() const
TString fInputAODGammaName
Input AOD gamma name.
TH2F ** fMixCOmega2
! MixC omega IVM(asy, pt, m)
Int_t fNBadChDistBin
Number of bad channel dist cut.
TH2F ** fRealOmega0
! Real omega IVM(asy, pt, m), with Asy_pi0<1
TH2F ** fMixBOmega1
! MixB omega IVM(asy, pt, m)
TH2F ** fMixCOmega1
! MixC omega IVM(asy, pt, m)
omega(782)->pi0+gamma->3gamma
Float_t GetHistoMassMax() const
Base class for CaloTrackCorr analysis algorithms.
Int_t fNVtxZBin
Number of vertex z cut.
int Int_t
Definition: External.C:63
virtual TClonesArray * GetInputAODBranch() const
virtual TClonesArray * GetAODMCParticles() const
virtual AliHistogramRanges * GetHistogramRanges()
Double_t * fCent
[fNVtxZBin]
virtual ~AliAnaOmegaToPi0Gamma()
Destructor.
float Float_t
Definition: External.C:68
TH2F ** fMixAOmega1
! MixA omega IVM(asy, pt, m)
const Double_t ptmax
virtual TClonesArray * GetAODBranch(const TString &aodBranchName) const
Recover ouput and input AOD pointers for each event in AliCaloTrackMaker.
virtual Int_t GetCentralityOpt() const
AliAnaOmegaToPi0Gamma()
Default constructor.
const Double_t ptmin
Bool_t IsBadRun(Int_t) const
Tests if this run bad according to private list.
virtual Double_t GetEventWeight() const
virtual TString GetAddedHistogramsStringToName() const
virtual Int_t GetNMaxEvMix() const
Number of bins in track multiplicity.
TH2F ** fMixBOmega0
! MixB omega IVM(asy, pt, m)
TH2F ** fRealOmega2
! Real omega IVM(asy, pt, m), with Asy_pi0<0.8
TH2F ** fMixAOmega2
! MixA omega IVM(asy, pt, m)
Int_t GetHistoPtBins() const
Double_t fPi0MassWindow
pi0 mass windows
TH2F ** fMixBOmega2
! MixB omega IVM(asy, pt, m)
void MakeAnalysisFillHistograms()
Do analysis, fill histograms.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Int_t fNCentBin
Number of centrality cut.
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
TClonesArray * fInputAODPi0
Input AOD pi0 array.
void Print(const Option_t *opt) const
Print some relevant parameters set in the analysis.
virtual AliCaloTrackReader * GetReader() const
TH1F * fhOmegaPriPt
! MC primary omega pt in 2pi and |y|<0.5
Double_t fPi0Mass
[fNRpBin]
Double_t fPi0OverOmegaPtCut
pi0 Pt over omega pt cut
TH2F * fhEtalon
! An etalon of 3D histograms
Double_t fEOverlapCluster
the pt when the two photons overlapped
Int_t nptbins
TH2F ** fMixAOmega0
! MixA omega IVM(asy, pt, m)
Double_t * fVtxZCut
Vertex z cut.
Int_t fNpid
Number of PID cut.
Double_t * fRp
[fNCentBin]
TList ** fEventsList
event list for mixing
void ReadHistograms(TList *outputList)