82addeb653ca4dc374cdc4ae67137e054967ea5c
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaOmegaToPi0Gamma.cxx
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 //_________________________________________________________________________
17 // class to extract omega(782)->pi0+gamma->3gamma
18 //  Mar. 22, 2011: Additional method, espeically for EMCAL. A high E cluster 
19 //  is assumpted as pi0 (two photons are overlapped) without unfolding
20 //
21 //-- Author: Renzhuo Wan (IOPP-Wuhan, China)
22 //_________________________________________________________________________
23
24 // --- ROOT system
25 class 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"
43 ClassImp(AliAnaOmegaToPi0Gamma)
44
45 //______________________________________________________________________________
46 AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma() : AliAnaCaloTrackCorrBaseClass(),
47 fInputAODPi0(0), fInputAODGammaName(""),
48 fEventsList(0x0),fNVtxZBin(0), fNCentBin(0), fNRpBin(0), fNBadChDistBin(0), fNpid(0),
49 fVtxZCut(0), fCent(0), fRp(0), 
50 fPi0Mass(0),fPi0MassWindow(0),fPi0OverOmegaPtCut(0),
51 fGammaOverOmegaPtCut(0), fEOverlapCluster(0),
52 fhEtalon(0),
53 fRealOmega0(0), fMixAOmega0(0),
54 fMixBOmega0(0), fMixCOmega0(0),
55 fRealOmega1(0), fMixAOmega1(0),
56 fMixBOmega1(0), fMixCOmega1(0),
57 fRealOmega2(0), fMixAOmega2(0),
58 fMixBOmega2(0), fMixCOmega2(0),
59 fhFakeOmega(0),
60 fhOmegaPriPt(0)
61 {
62  //Default Ctor
63  InitParameters();
64 }
65
66 //______________________________________________________________________________
67 AliAnaOmegaToPi0Gamma::~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 //______________________________________________________________________________
96 void 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;                   
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 //______________________________________________________________________________
117 TList * 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   //
129   Int_t nptbins   = GetHistogramRanges()->GetHistoPtBins();
130   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();
131   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();
132   
133   Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
134   Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
135   Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
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 //______________________________________________________________________________
296 void 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() ) ;
300   AliAnaCaloTrackCorrBaseClass::Print(" ");
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);
309
310
311 //______________________________________________________________________________
312 void 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();
347       if(!mcparticles0 )
348       {
349         if(GetDebug() > 0) printf("AliAnaAcceptance::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
350       }
351       else
352       {
353         for(Int_t i=0;i<mcparticles0->GetEntries();i++)
354         {
355           aodprimary =(AliAODMCParticle*)mcparticles0->At(i);
356           pdg = aodprimary->GetPdgCode() ;
357           eta=aodprimary->Eta();
358           pt=aodprimary->Pt();
359           if(TMath::Abs(eta)<0.5)
360           {
361             if(pdg==223) fhOmegaPriPt->Fill(pt);
362           }
363           
364         }
365       }//mcparticles0 exists
366     }//AOD MC Particles
367   }// is data and MC
368   
369   
370   //process event from AOD brach 
371   //extract pi0, eta and omega analysis
372   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
373   if(IsBadRun(iRun)) return ;   
374   
375   //vertex z
376   Double_t vert[]={0,0,0} ;
377   GetVertex(vert);
378   Int_t curEventBin =0;
379   
380   Int_t ivtxzbin=(Int_t)TMath::Abs(vert[2])/10;
381   if(ivtxzbin>=fNVtxZBin)return;
382   
383   //centrality
384   Int_t currentCentrality = GetEventCentrality();
385   if(currentCentrality == -1) return;
386   Int_t optCent = GetReader()->GetCentralityOpt();
387   Int_t icentbin=currentCentrality/(optCent/fNCentBin) ; //GetEventCentrality();
388  
389   printf("-------------- %d  %d  %d  ",currentCentrality, optCent, icentbin);
390   //reaction plane
391   Int_t irpbin=0;
392   
393   if(ivtxzbin==-1) return; 
394   curEventBin = ivtxzbin*fNCentBin*fNRpBin + icentbin*fNRpBin + irpbin;
395   TClonesArray *aodGamma = (TClonesArray*) GetAODBranch(fInputAODGammaName); //photon array
396   //  TClonesArray *aodGamma = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaName); //photon array
397   Int_t nphotons =0;
398   if(aodGamma) nphotons= aodGamma->GetEntries(); 
399   else return;
400   
401   fInputAODPi0 = (TClonesArray*)GetInputAODBranch();  //pi0 array
402   Int_t npi0s = 0;
403   if(fInputAODPi0) npi0s= fInputAODPi0 ->GetEntries();
404   else return;
405   
406   if(nphotons<3 || npi0s<1)return; //for pi0, eta and omega->pi0+gamma->3gamma reconstruction
407   
408   //reconstruction of omega(782)->pi0+gamma->3gamma
409   //loop for pi0 and photon
410   if(GetDebug() > 0) printf("omega->pi0+gamma->3gamma invariant mass analysis ! This event have %d photons and %d pi0 \n", nphotons, npi0s);
411   for(Int_t i=0;i<npi0s;i++){
412     AliAODPWG4Particle * pi0 = (AliAODPWG4Particle*) (fInputAODPi0->At(i)) ; //pi0
413     TLorentzVector vpi0(pi0->Px(),pi0->Py(),pi0->Pz(),pi0->E());
414     Int_t lab1=pi0->GetCaloLabel(0);  // photon1 from pi0 decay
415     Int_t lab2=pi0->GetCaloLabel(1);  // photon2 from pi0 decay
416     //for omega->pi0+gamma, it needs at least three photons per event
417     //Get the two decay photons from pi0
418     AliAODPWG4Particle * photon1 =0;
419     AliAODPWG4Particle * photon2 =0;
420     for(Int_t d1=0;d1<nphotons;d1++){
421       for(Int_t d2=0;d2<nphotons;d2++){
422         AliAODPWG4Particle * dp1 = (AliAODPWG4Particle*) (aodGamma->At(d1));
423         AliAODPWG4Particle * dp2 = (AliAODPWG4Particle*) (aodGamma->At(d2));
424         Int_t dlab1=dp1->GetCaloLabel(0);
425         Int_t dlab2=dp2->GetCaloLabel(0);
426         if(dlab1==lab1 && dlab2==lab2){
427           photon1=dp1;
428           photon2=dp2;
429         }
430         else continue;
431       }
432     }
433     //caculate the asy and dist of the two photon from pi0 decay
434     TLorentzVector dph1(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());
435     TLorentzVector dph2(photon2->Px(),photon2->Py(),photon2->Pz(),photon2->E());
436     
437     Double_t pi0asy= TMath::Abs(dph1.E()-dph2.E())/(dph1.E()+dph2.E());
438     //    Double_t phi1=dph1.Phi();
439     //    Double_t phi2=dph2.Phi();
440     //    Double_t eta1=dph1.Eta();
441     //    Double_t eta2=dph2.Eta();
442     //    Double_t pi0dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));
443     
444     if(pi0->GetIdentifiedParticleType()==111  && nphotons>2 && npi0s
445        && TMath::Abs(vpi0.M()-fPi0Mass)<fPi0MassWindow) { //pi0 candidates
446       
447       //avoid the double counting
448       Int_t * dc1= new Int_t[nphotons];
449       Int_t * dc2= new Int_t[nphotons];
450       Int_t index1=0;
451       Int_t index2=0;
452       for(Int_t k=0;k<i;k++){
453         AliAODPWG4Particle * p3=(AliAODPWG4Particle*)(fInputAODPi0->At(k));
454         Int_t lab4=p3->GetCaloLabel(0);
455         Int_t lab5=p3->GetCaloLabel(1);
456         if(lab1==lab4){ dc1[index1]=lab5;  index1++;  }
457         if(lab2==lab5){ dc2[index2]=lab4;  index2++;  }
458       }
459       
460       
461       //loop the pi0 with third gamma
462       for(Int_t j=0;j<nphotons;j++){
463         AliAODPWG4Particle *photon3 = (AliAODPWG4Particle*) (aodGamma->At(j));
464         TLorentzVector dph3(photon3->Px(),photon3->Py(),photon3->Pz(),photon3->E());
465         Int_t lab3=photon3->GetCaloLabel(0);
466         Double_t pi0gammapt=(vpi0+dph3).Pt();
467         Double_t pi0gammamass=(vpi0+dph3).M();
468         Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt; 
469         Double_t gammaOverOmegaPtRatio= dph3.Pt()/pi0gammapt;
470         
471         //pi0, gamma pt cut             
472         if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut || 
473            gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
474         
475         for(Int_t l=0;l<index1;l++) if(lab3==dc1[l]) lab3=-1;
476         for(Int_t l=0;l<index2;l++) if(lab3==dc2[l]) lab3=-1;
477         
478         if(lab3>0 && lab3!=lab1 && lab3!=lab2){
479                 for(Int_t ipid=0;ipid<fNpid;ipid++){
480             for(Int_t idist=0;idist<fNBadChDistBin;idist++){
481               Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
482               if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
483                  photon2->IsPIDOK(ipid,AliCaloPID::kPhoton) && 
484                  photon3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
485                  photon1->DistToBad()>=idist &&
486                  photon2->DistToBad()>=idist &&
487                  photon3->DistToBad()>=idist ){
488                 //fill the histograms
489                 if(GetDebug() > 2) printf("Real: index  %d  pt  %2.3f  mass   %2.3f \n", index, pi0gammapt, pi0gammamass);
490                 fRealOmega0[index]->Fill(pi0gammapt,pi0gammamass);
491                 if(pi0asy<0.7) fRealOmega1[index]->Fill(pi0gammapt,pi0gammamass);
492                 if(pi0asy<0.8) fRealOmega2[index]->Fill(pi0gammapt,pi0gammamass);
493               }
494             }
495                 }
496         }
497       } 
498       delete []dc1;
499       delete []dc2;
500       if(GetDebug() > 0) printf("MixA: (r1_event1+r2_event1)+r3_event2 \n");
501       //-------------------------
502       //background analysis
503       //three background
504       // --A   (r1_event1+r2_event1)+r3_event2
505       Int_t nMixed = fEventsList[curEventBin]->GetSize();
506       for(Int_t im=0;im<nMixed;im++){
507         TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(im));
508         for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
509           AliAODPWG4Particle *mix1ph = (AliAODPWG4Particle*) (ev2->At(mix1));     
510           TLorentzVector vmixph(mix1ph->Px(),mix1ph->Py(),mix1ph->Pz(),mix1ph->E());
511           Double_t pi0gammapt=(vpi0+vmixph).Pt();
512           Double_t pi0gammamass=(vpi0+vmixph).M();
513           Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
514           Double_t gammaOverOmegaPtRatio= vmixph.Pt()/pi0gammapt;
515           
516           //pi0, gamma pt cut             
517           if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut || 
518              gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
519           
520                 for(Int_t ipid=0;ipid<fNpid;ipid++){
521             for(Int_t idist=0;idist<fNBadChDistBin;idist++){
522               Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
523               if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
524                  photon2->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
525                  mix1ph->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
526                  photon1->DistToBad()>=idist &&
527                  photon2->DistToBad()>=idist &&
528                  mix1ph->DistToBad()>=idist ){
529                 if(GetDebug() > 2) printf("MixA: index  %d   pt  %2.3f  mass   %2.3f \n",index, pi0gammapt, pi0gammamass);
530                 //fill the histograms
531                 fMixAOmega0[index]->Fill(pi0gammapt,pi0gammamass);
532                 if(pi0asy<0.7)fMixAOmega1[index]->Fill(pi0gammapt,pi0gammamass);
533                 if(pi0asy<0.8)fMixAOmega2[index]->Fill(pi0gammapt,pi0gammamass);
534                 //printf("mix A  %d  %2.2f \n", index, pi0gammamass);
535                 
536               }
537             }
538           }
539         }
540       }
541     }
542   }
543   
544   //
545   // --B   (r1_event1+r2_event2)+r3_event2
546   //
547   if(GetDebug() >0)printf("MixB:  (r1_event1+r2_event2)+r3_event2 \n");
548   for(Int_t i=0;i<nphotons;i++){
549     AliAODPWG4Particle *ph1 = (AliAODPWG4Particle*) (aodGamma->At(i)); 
550     TLorentzVector vph1(ph1->Px(),ph1->Py(),ph1->Pz(),ph1->E());
551     //interrupt here...................
552     //especially for EMCAL
553     //we suppose the high pt clusters are overlapped pi0   
554
555     for(Int_t j=i+1;j<nphotons;j++){
556         AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (aodGamma->At(j));
557         TLorentzVector fakePi0, fakeOmega, vph;
558
559         if(ph1->E() > fEOverlapCluster && ph1->E() > ph2->E()) {
560            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));
561            vph.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
562         }
563         else if(ph2->E() > fEOverlapCluster && ph2->E() > ph1->E()) {
564            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));
565            vph.SetPxPyPzE(ph1->Px(), ph1->Py(),ph1->Pz(),ph1->E());
566         }
567         else continue;
568
569         fakeOmega=fakePi0+vph;
570         for(Int_t ii=0;ii<fNCentBin;ii++){ 
571            fhFakeOmega[icentbin]->Fill(fakeOmega.Pt(), fakeOmega.M());
572         }
573     }//j
574
575     //continue ................
576     Int_t nMixed = fEventsList[curEventBin]->GetSize();
577     for(Int_t ie=0;ie<nMixed;ie++){
578       TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(ie));
579       for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
580         AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (ev2->At(mix1));
581         TLorentzVector vph2(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
582         Double_t pi0asy = TMath::Abs(vph1.E()-vph2.E())/(vph1.E()+vph2.E());         
583         Double_t pi0mass=(vph1+vph2).M();
584         
585         if(TMath::Abs(pi0mass-fPi0Mass)<fPi0MassWindow){//for pi0 selection
586           for(Int_t mix2=(mix1+1);mix2<ev2->GetEntries();mix2++){
587             AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev2->At(mix2));
588             TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
589             
590             Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
591             Double_t pi0gammamass=(vph1+vph2+vph3).M(); 
592             Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
593             Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
594             //pi0, gamma pt cut             
595             if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
596                gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
597             
598             for(Int_t ipid=0;ipid<fNpid;ipid++){
599               for(Int_t idist=0;idist<fNBadChDistBin;idist++){
600                 Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
601                 if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
602                                ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
603                    ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
604                    ph1->DistToBad()>=idist &&
605                    ph2->DistToBad()>=idist &&
606                    ph3->DistToBad()>=idist ){
607                   if(GetDebug() > 2) printf("MixB: index  %d   pt  %2.3f  mass   %2.3f \n", index, pi0gammapt, pi0gammamass);
608                   //fill histograms
609                   fMixBOmega0[index]->Fill(pi0gammapt,pi0gammamass);
610                   if(pi0asy<0.7) fMixBOmega1[index]->Fill(pi0gammapt,pi0gammamass);
611                   if(pi0asy<0.8) fMixBOmega2[index]->Fill(pi0gammapt,pi0gammamass);
612                   //printf("mix B  %d  %2.2f \n", index, pi0gammamass);
613                 }
614               }             
615             }
616           }
617           
618           //
619                 // --C   (r1_event1+r2_event2)+r3_event3
620           //
621           if(GetDebug() >0)printf("MixC: (r1_event1+r2_event2)+r3_event3\n");
622           for(Int_t je=(ie+1);je<nMixed;je++){
623             TClonesArray* ev3= (TClonesArray*) (fEventsList[curEventBin]->At(je));
624             for(Int_t mix3=0;mix3<ev3->GetEntries();mix3++){
625               AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev3->At(mix3));
626               TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
627               
628               Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
629               Double_t pi0gammamass=(vph1+vph2+vph3).M();
630               Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
631               Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
632               //pi0, gamma pt cut             
633               if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
634                  gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
635               
636               for(Int_t ipid=0;ipid<fNpid;ipid++){
637                 for(Int_t idist=0;idist<fNBadChDistBin;idist++){
638                   Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
639                   if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
640                      ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
641                      ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
642                      ph1->DistToBad()>=idist &&
643                      ph2->DistToBad()>=idist &&
644                      ph3->DistToBad()>=idist ){
645                     if(GetDebug() > 2) printf("MixC: index  %d  pt  %2.3f  mass   %2.3f \n", index, pi0gammapt, pi0gammamass);
646                     //fill histograms
647                     fMixCOmega0[index]->Fill(pi0gammapt,pi0gammamass);
648                     if(pi0asy<0.7) fMixCOmega1[index]->Fill(pi0gammapt,pi0gammamass);
649                     if(pi0asy<0.8) fMixCOmega2[index]->Fill(pi0gammapt,pi0gammamass);
650                     //printf("mix C  %d  %2.2f \n", index, pi0gammamass);
651                   }
652                 }
653               }
654             }
655           }
656         } //for pi0 selecton            
657       }
658     }
659   }
660   
661   
662   //event buffer 
663   TClonesArray *currentEvent = new TClonesArray(*aodGamma);
664   if(currentEvent->GetEntriesFast()>0){
665     fEventsList[curEventBin]->AddFirst(currentEvent) ;
666     currentEvent=0 ; 
667     if(fEventsList[curEventBin]->GetSize()>=GetNMaxEvMix()) {
668       TClonesArray * tmp = (TClonesArray*) (fEventsList[curEventBin]->Last()) ;
669       fEventsList[curEventBin]->RemoveLast() ;
670       delete tmp ;
671     }
672   }
673   else{ 
674     delete currentEvent ;
675     currentEvent=0 ;
676   }
677   
678 }
679
680 //______________________________________________________________________________
681 void AliAnaOmegaToPi0Gamma::ReadHistograms(TList * outputList)
682 {
683  //read the histograms 
684  //for the finalization of the terminate analysis
685
686  Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"RealToPi0Gamma_Vz0C0Rp0Pid0Dist0"));
687
688   Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
689
690  if(!fRealOmega0) fRealOmega0 =new TH2F*[ndim];
691  if(!fMixAOmega0) fMixAOmega0 =new TH2F*[ndim];
692  if(!fMixBOmega0) fMixBOmega0 =new TH2F*[ndim];
693  if(!fMixCOmega0) fMixCOmega0 =new TH2F*[ndim];
694
695  if(!fRealOmega1) fRealOmega1 =new TH2F*[ndim];
696  if(!fMixAOmega1) fMixAOmega1 =new TH2F*[ndim];
697  if(!fMixBOmega1) fMixBOmega1 =new TH2F*[ndim];
698  if(!fMixCOmega1) fMixCOmega1 =new TH2F*[ndim];
699
700  if(!fRealOmega2) fRealOmega2 =new TH2F*[ndim];
701  if(!fMixAOmega2) fMixAOmega2 =new TH2F*[ndim];
702  if(!fMixBOmega2) fMixBOmega2 =new TH2F*[ndim];
703  if(!fMixCOmega2) fMixCOmega2 =new TH2F*[ndim];
704
705   for(Int_t i=0;i<fNVtxZBin;i++){
706      for(Int_t j=0;j<fNCentBin;j++){
707          for(Int_t k=0;k<fNRpBin;k++){ //at event level
708              Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
709              for(Int_t ipid=0;ipid<fNpid;ipid++){ 
710                 for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle
711                     Int_t ind=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
712                     fRealOmega0[ind]= (TH2F*) outputList->At(index++);
713                     fMixAOmega0[ind]= (TH2F*) outputList->At(index++);
714                     fMixBOmega0[ind]= (TH2F*) outputList->At(index++);
715                     fMixCOmega0[ind]= (TH2F*) outputList->At(index++);
716
717                     fRealOmega1[ind]= (TH2F*) outputList->At(index++);
718                     fMixAOmega1[ind]= (TH2F*) outputList->At(index++);
719                     fMixBOmega1[ind]= (TH2F*) outputList->At(index++);
720                     fMixCOmega1[ind]= (TH2F*) outputList->At(index++);
721
722                     fRealOmega2[ind]= (TH2F*) outputList->At(index++);
723                     fMixAOmega2[ind]= (TH2F*) outputList->At(index++);
724                     fMixBOmega2[ind]= (TH2F*) outputList->At(index++);
725                     fMixCOmega2[ind]= (TH2F*) outputList->At(index++);
726                     
727                  
728                 }
729               }
730           }
731       }
732   }
733   
734   if(IsDataMC()){
735      fhOmegaPriPt  = (TH1F*)  outputList->At(index++);
736   }
737
738 }
739
740 //______________________________________________________________________________
741 void AliAnaOmegaToPi0Gamma::Terminate(TList * outputList) 
742 {
743 // //Do some calculations and plots from the final histograms.
744   if(GetDebug() >= 0) printf("AliAnaOmegaToPi0Gamma::Terminate() \n");
745   ReadHistograms(outputList);
746   const Int_t buffersize = 255;
747   char cvs1[buffersize];  
748   snprintf(cvs1, buffersize, "Neutral_%s_IVM",fInputAODGammaName.Data());
749
750   TCanvas * cvsIVM = new TCanvas(cvs1, cvs1, 400, 10, 600, 700) ;
751   cvsIVM->Divide(2, 2);
752
753   cvsIVM->cd(1);
754   char dec[buffersize];
755   snprintf(dec,buffersize,"h2Real_%s",fInputAODGammaName.Data());
756   TH2F * h2Real= (TH2F*)fRealOmega0[0]->Clone(dec);
757   h2Real->GetXaxis()->SetRangeUser(4,6);
758   TH1F * hRealOmega = (TH1F*) h2Real->ProjectionY();
759   hRealOmega->SetTitle("RealPi0Gamma 4<pt<6");
760   hRealOmega->SetLineColor(2);
761   hRealOmega->Draw();
762
763   cvsIVM->cd(2);
764   snprintf(dec,buffersize,"hMixA_%s",fInputAODGammaName.Data());
765   TH2F *h2MixA= (TH2F*)fMixAOmega0[0]->Clone(dec);
766   h2MixA->GetXaxis()->SetRangeUser(4,6);
767   TH1F * hMixAOmega = (TH1F*) h2MixA->ProjectionY();
768   hMixAOmega->SetTitle("MixA 4<pt<6");
769   hMixAOmega->SetLineColor(2);
770   hMixAOmega->Draw();
771
772   cvsIVM->cd(3);
773   snprintf(dec,buffersize,"hMixB_%s",fInputAODGammaName.Data());
774   TH2F * h2MixB= (TH2F*)fMixBOmega0[0]->Clone(dec);
775   h2MixB->GetXaxis()->SetRangeUser(4,6);
776   TH1F * hMixBOmega = (TH1F*) h2MixB->ProjectionY();
777   hMixBOmega->SetTitle("MixB 4<pt<6");
778   hMixBOmega->SetLineColor(2);
779   hMixBOmega->Draw();
780
781   cvsIVM->cd(4);
782   snprintf(dec,buffersize,"hMixC_%s",fInputAODGammaName.Data());
783   TH2F *h2MixC= (TH2F*)fMixCOmega0[0]->Clone(dec);
784   h2MixC->GetXaxis()->SetRangeUser(4,6);
785   TH1F * hMixCOmega = (TH1F*) h2MixC->ProjectionY();
786   hMixCOmega->SetTitle("MixC 4<pt<6");
787   hMixCOmega->SetLineColor(2);
788   hMixCOmega->Draw();
789
790   char eps[buffersize];
791   snprintf(eps,buffersize,"CVS_%s_IVM.eps",fInputAODGammaName.Data());
792   cvsIVM->Print(eps);
793   cvsIVM->Modified();
794  
795 }