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