Jet and Particle identification tasks moved to different directories
[u/mrichter/AliRoot.git] / PWG4 / PartCorr / AliAnaPi0.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 /* $Id: $ */
16
17 //_________________________________________________________________________
18 // Class to collect two-photon invariant mass distributions for
19 // extractin raw pi0 yield.
20 //
21 //-- Author: Dmitri Peressounko (RRC "KI") 
22 //_________________________________________________________________________
23
24
25 // --- ROOT system ---
26 #include "TClonesArray.h"
27 #include "TRefArray.h"
28 #include "TH3.h"
29
30 //---- AliRoot system ----
31 #include "AliAnaPi0.h"
32 #include "AliLog.h"
33 #include "AliCaloPhoton.h"
34 #include "AliESDEvent.h"
35 #include "AliAODEvent.h"
36 #include "AliESDVertex.h"
37 #include "AliAODVertex.h"
38 #include "AliESDCaloCluster.h"
39 #include "AliAODCaloCluster.h"
40
41 ClassImp(AliAnaPi0)
42   
43 //____________________________________________________________________________
44   AliAnaPi0::AliAnaPi0() : AliAnalysisTaskSE(),
45   fNCentrBin(1),fNZvertBin(1),fNrpBin(1),fNPID(2),fNmaxMixEv(10),
46   fCurCentrBin(0),fCurZvertBin(0),fCurRPBin(0),fPtMin(0.),fZvtxCut(40.),
47   fMinDist(2.),fMinDist2(4.),fMinDist3(5.),fDispCut(1.5),fTOFCut(5.e-9),fPhotPID(0.6),
48   fEventsList(0x0),fCurrentEvent(0x0),fOutputList(0x0),fhEtalon(0x0),
49   fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0) 
50 {
51   //Default Ctor
52
53 }
54 //____________________________________________________________________________
55   AliAnaPi0::AliAnaPi0(const char *name): AliAnalysisTaskSE(name),
56   fNCentrBin(1),fNZvertBin(1),fNrpBin(1),fNPID(9),fNmaxMixEv(10),
57   fCurCentrBin(0),fCurZvertBin(0),fCurRPBin(0),fPtMin(0.),fZvtxCut(40.),
58   fMinDist(2.),fMinDist2(4.),fMinDist3(5.),fDispCut(1.5),fTOFCut(5.e-9),fPhotPID(0.6),
59   fEventsList(0x0),fCurrentEvent(0x0),fOutputList(0x0),fhEtalon(0x0),
60   fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0) 
61 {
62   //Ctor
63   DefineOutput(1,TList::Class());
64 }
65
66 //____________________________________________________________________________
67 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnalysisTaskSE(ex),  
68   fNCentrBin(1),fNZvertBin(1),fNrpBin(1),fNPID(9),fNmaxMixEv(10),
69   fCurCentrBin(0),fCurZvertBin(0),fCurRPBin(0),fPtMin(0.),fZvtxCut(40.),
70   fMinDist(2.),fMinDist2(4.),fMinDist3(5.),fDispCut(1.5),fTOFCut(5.e-9),fPhotPID(0.6),
71   fEventsList(0x0),fCurrentEvent(0x0),fOutputList(0x0),fhEtalon(0x0),
72   fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0) 
73 {
74   // cpy ctor
75   //Do not need it
76 }
77 //_________________________________________________________________________
78 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
79 {
80   // assignment operator
81
82   if(this == &ex)return *this;
83   ((AliAnalysisTaskSE *)this)->operator=(ex);
84  
85   return *this;
86
87 }
88 //____________________________________________________________________________
89 AliAnaPi0::~AliAnaPi0() {
90   // Remove event containers
91   if(fEventsList){
92     for(Int_t ic=0; ic<fNCentrBin; ic++){
93       for(Int_t iz=0; iz<fNZvertBin; iz++){
94         for(Int_t irp=0; irp<fNrpBin; irp++){
95           fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
96           delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
97         }
98       }
99     }
100     delete[] fEventsList; 
101     fEventsList=0 ;
102   }
103
104   if(fhEtalon){
105     delete fhEtalon ;
106     fhEtalon = 0; 
107   }
108   if(fhRe1) delete[] fhRe1 ; //Do not delete histograms!
109   if(fhRe2) delete[] fhRe2 ; //Do not delete histograms!
110   if(fhRe3) delete[] fhRe3 ; //Do not delete histograms!
111   if(fhMi1) delete[] fhMi1 ; //Do not delete histograms!
112   if(fhMi2) delete[] fhMi2 ; //Do not delete histograms!
113   if(fhMi3) delete[] fhMi3 ; //Do not delete histograms!
114
115 }
116 //________________________________________________________________________
117 void  AliAnaPi0::UserCreateOutputObjects()
118 {  
119   // Create histograms to be saved in output file and 
120   // store them in fOutputContainer
121
122   AliDebug(1,"Init inv. mass histograms");
123   OpenFile(1) ; 
124
125   //create event containers
126   fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
127   for(Int_t ic=0; ic<fNCentrBin; ic++){
128     for(Int_t iz=0; iz<fNZvertBin; iz++){
129       for(Int_t irp=0; irp<fNrpBin; irp++){
130         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
131       }
132     }
133   }
134
135   fOutputList = new TList() ;
136   fOutputList->SetName(GetName()) ;
137  
138   fhRe1=new TH3D*[fNCentrBin*fNPID] ;
139   fhRe2=new TH3D*[fNCentrBin*fNPID] ;
140   fhRe3=new TH3D*[fNCentrBin*fNPID] ;
141   fhMi1=new TH3D*[fNCentrBin*fNPID] ;
142   fhMi2=new TH3D*[fNCentrBin*fNPID] ;
143   fhMi3=new TH3D*[fNCentrBin*fNPID] ;
144
145   char key[255] ;
146   char title[255] ;
147   for(Int_t ic=0; ic<fNCentrBin; ic++){
148     for(Int_t ipid=0; ipid<fNPID; ipid++){
149       //Distance to bad module 1
150       sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
151       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
152       fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
153       fhRe1[ic*fNPID+ipid]->SetName(key) ;
154       fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
155       fOutputList->Add(fhRe1[ic*fNPID+ipid]) ;
156
157       sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
158       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
159       fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
160       fhMi1[ic*fNPID+ipid]->SetName(key) ;
161       fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
162       fOutputList->Add(fhMi1[ic*fNPID+ipid]) ;
163
164       //Distance to bad module 2
165       sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
166       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
167       fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
168       fhRe2[ic*fNPID+ipid]->SetName(key) ;
169       fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
170       fOutputList->Add(fhRe2[ic*fNPID+ipid]) ;
171
172       sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
173       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
174       fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
175       fhMi2[ic*fNPID+ipid]->SetName(key) ;
176       fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
177       fOutputList->Add(fhMi2[ic*fNPID+ipid]) ;
178
179       //Distance to bad module 3
180       sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
181       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
182       fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
183       fhRe3[ic*fNPID+ipid]->SetName(key) ; 
184       fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
185       fOutputList->Add(fhRe3[ic*fNPID+ipid]) ;
186
187       sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
188       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
189       fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
190       fhMi3[ic*fNPID+ipid]->SetName(key) ;
191       fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
192       fOutputList->Add(fhMi3[ic*fNPID+ipid]) ;
193     }
194   }
195
196   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
197                     fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
198   fOutputList->Add(fhEvents) ;
199
200   //Save parameters used for analysis
201   TString parList ; //this will be list of parameters used for this analysis.
202   char onePar[255] ;
203   sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
204   parList+=onePar ;
205   sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
206   parList+=onePar ;
207   sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
208   parList+=onePar ;
209   sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
210   parList+=onePar ;
211   sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
212   parList+=onePar ;
213   sprintf(onePar,"Cuts: \n") ;
214   parList+=onePar ;
215   sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
216   parList+=onePar ;
217   sprintf(onePar,"Minimal P_t: %f \n", fPtMin) ;
218   parList+=onePar ;
219   sprintf(onePar,"fMinDist =%f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
220   parList+=onePar ;
221   sprintf(onePar,"fMinDist2=%f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
222   parList+=onePar ;
223   sprintf(onePar,"fMinDist3=%f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
224   parList+=onePar ;
225   sprintf(onePar,"fDispCut =%f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
226   parList+=onePar ;
227   sprintf(onePar,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
228   parList+=onePar ;
229   sprintf(onePar,"fPhotPID =%f (limit for Baesian PID for photon) \n",fPhotPID) ;
230   parList+=onePar ;
231  
232   TObjString *oString= new TObjString(parList) ;
233   fOutputList->Add(oString);
234  
235 }
236 //__________________________________________________
237 void AliAnaPi0::InitParameters()
238
239   //Initialize the parameters of the analysis.
240   
241 }
242  //__________________________________________________
243 void AliAnaPi0::Init()
244
245   //Make here all memory allocations
246   //create etalon histo for all later histograms
247   if(!fhEtalon){                                                   //  p_T      alpha   d m_gg    
248     fhEtalon = new TH3D("hEtalon","Histo with binning parameters",20,0.,10.,10,0.,1.,200,0.,1.) ; 
249     fhEtalon->SetXTitle("P_{T} (GeV)") ;
250     fhEtalon->SetYTitle("#alpha") ;
251     fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
252   }
253   
254 }
255 //__________________________________________________________________
256 void AliAnaPi0::Print(const Option_t * /*opt*/) const
257 {
258   //Print some relevant parameters set for the analysis
259   printf("Class AliAnaPi0 for gamma-gamma inv.mass construction \n") ;
260   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
261   printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
262   printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
263   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
264   printf("Number of different PID used:  %d \n",fNPID) ;
265   printf("Cuts: \n") ;
266   printf("Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
267   printf("Minimal P_t: %f \n", fPtMin) ;
268   printf("fMinDist =%f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
269   printf("fMinDist2=%f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
270   printf("fMinDist3=%f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
271   printf("fDispCut =%f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
272   printf("fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
273   printf("fPhotPID =%f (limit for Baesian PID for photon) \n",fPhotPID) ;
274   printf("------------------------------------------------------\n") ;
275   
276
277 //__________________________________________________________________
278 void AliAnaPi0::UserExec(Option_t *)
279 {
280   //Process one event and extract photons 
281   //impose cuts on bad modules and bad runs
282   //fill internal storage for subsequent histo filling
283
284   AliVEvent* event = InputEvent();
285  
286   //Apply some cuts on event: vertex position and centrality range  
287   Int_t iRun=event->GetRunNumber() ;
288   if(IsBadRun(iRun))
289     return ;
290
291   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
292   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event) ;
293   if(esd){
294      if(!FillFromESD(esd))
295       return  ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
296    }
297    else if(aod){
298      if(!FillFromAOD(aod))
299        return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
300    }
301    else{
302      printf("Input neither ESD nor AOD, do nothing \n") ;
303      return ;
304    }
305
306
307   fhEvents->Fill(fCurCentrBin+0.5,fCurZvertBin+0.5,fCurRPBin+0.5) ;
308
309   if(fCurrentEvent->GetEntriesFast()>0){
310     //Reduce size for storing
311     fCurrentEvent->Expand(fCurrentEvent->GetEntriesFast()) ;
312     FillHistograms() ;
313   }
314
315   PostData(1, fOutputList);
316  
317
318 //__________________________________________________________________
319 Bool_t AliAnaPi0::FillFromESD(AliESDEvent * esd){
320   //Fill photon list from ESD applying 
321   //some cut should be applyed
322
323   //Impose cut on vertex and calculate Zvertex bin
324   const AliESDVertex *esdV = esd->GetVertex() ;
325   esdV->GetXYZ(fVert);
326   if(fVert[2]<-fZvtxCut || fVert[2]> fZvtxCut)
327     return kFALSE ;
328   fCurZvertBin=(Int_t)(0.5*fNZvertBin*(fVert[2]+fZvtxCut)/fZvtxCut) ;
329     
330   //Get Centrality and calculate centrality bin
331   //Does not exist in ESD yet???????
332   fCurCentrBin=0 ;
333  
334   //Get Reaction Plain position and calculate RP bin
335   //does not exist in ESD yet????
336   fCurRPBin=0 ;
337  
338   //************************  PHOS *************************************
339   TRefArray * caloClustersArr  = new TRefArray();
340   esd->GetPHOSClusters(caloClustersArr);
341  
342   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
343   //if fCurrentEvent!=0 it was not filled in previous event, still empty and can be used
344   if(!fCurrentEvent) 
345     fCurrentEvent = new TClonesArray("AliCaloPhoton",kNumberOfPhosClusters) ;
346   Int_t inList=0 ;
347  
348   // loop over the PHOS Cluster
349   for(Int_t i = 0 ; i < kNumberOfPhosClusters ; i++) {
350     AliESDCaloCluster * calo = (AliESDCaloCluster *) caloClustersArr->At(i) ;
351  
352     //Make some tests
353     if(calo->E()<fPtMin)
354       continue ;
355
356     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad in cm
357     if(distBad<0.)distBad=9999. ; //workout strange convension dist = -1. ;
358     if(distBad<fMinDist) //In bad channel (cristall size 2.2x2.2 cm)
359       continue ;
360
361     new((*fCurrentEvent)[inList])AliCaloPhoton() ;
362     AliCaloPhoton * ph = static_cast<AliCaloPhoton*>(fCurrentEvent->At(inList)) ;
363     inList++ ;
364
365     //Set Momentum 
366     calo->GetMomentum(*ph,fVert);
367
368     //Dispersion/lambdas
369     Double_t disp=calo->GetClusterDisp()  ;
370 //    Double_t m20=calo->GetM20() ;
371 //    Double_t m02=calo->GetM02() ; 
372     ph->SetDispBit(disp<fDispCut) ;  
373
374     //TOF
375     Double_t tof=calo->GetTOF()  ;
376     ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ; 
377  
378     //Charged veto
379 //    Double_t cpvR=calo->GetEmcCpvDistance() ; 
380     Int_t itr=calo->GetNTracksMatched();  //number of track matched
381     ph->SetCPVBit(itr>0) ;  //Temporary cut, should we evaluate distance?
382
383     //Overall PID
384     Double_t *pid=calo->GetPid();
385     ph->SetPCAPID(pid[AliPID::kPhoton]>fPhotPID) ;
386       
387     //Set Distance to Bad channel
388     if(distBad>fMinDist3)
389       ph->SetDistToBad(2) ;
390     else if(distBad>fMinDist2)
391         ph->SetDistToBad(1) ; 
392     else
393        ph->SetDistToBad(0) ;
394
395   }//loop
396   return kTRUE ;
397 }
398 //__________________________________________________________________
399 Bool_t AliAnaPi0::FillFromAOD(AliAODEvent * aod){
400   //Fill photon list from AOD applying 
401   //some cuts
402
403   //Impose cut on vertex and calculate Zvertex bin
404   const AliAODVertex *aodV = aod->GetPrimaryVertex() ;
405   fVert[0]=aodV->GetX();
406   fVert[1]=aodV->GetY();
407   fVert[2]=aodV->GetZ();
408   if(fVert[2]<-fZvtxCut || fVert[2]> fZvtxCut)
409     return kFALSE ;
410   fCurZvertBin=(Int_t)(0.5*fNZvertBin*(fVert[2]+fZvtxCut)/fZvtxCut) ;
411  
412   //Get Centrality and calculate centrality bin
413   //Does not exist in ESD yet???????
414   fCurCentrBin=0 ;
415  
416   //Get Reaction Plain position and calculate RP bin
417   //does not exist in ESD yet????
418   fCurRPBin=0 ;
419  
420   //************************  PHOS *************************************
421   const Int_t kNumberOfPhosClusters   = aod->GetNCaloClusters() ;
422   if(!fCurrentEvent)
423     fCurrentEvent = new TClonesArray("AliCaloPhoton",kNumberOfPhosClusters) ;
424   Int_t inList=0 ;
425  
426   // loop over the PHOS Cluster
427   for(Int_t i = 0 ; i < kNumberOfPhosClusters ; i++) {
428     AliAODCaloCluster * calo = aod->GetCaloCluster(i);
429       
430     //Make some tests
431     if(calo->E()<fPtMin)
432       continue ;
433
434     Double_t distBad=calo->GetDistToBadChannel() ;
435     if(distBad<fMinDist) //In bad channel
436       continue ;
437
438     new((*fCurrentEvent)[inList])AliCaloPhoton() ;
439     AliCaloPhoton * ph = static_cast<AliCaloPhoton*>(fCurrentEvent->At(inList)) ;
440     inList++ ;
441
442     //Set Momentum 
443     calo->GetMomentum(*ph,fVert);
444
445     //Dispersion/lambdas
446     Double_t disp=calo->GetDispersion()  ;
447 //    Double_t m20=calo->GetM20() ;
448 //    Double_t m02=calo->GetM02() ; 
449     ph->SetDispBit(disp<fDispCut) ; 
450
451     //TOF
452     Double_t tof=calo->GetTOF()  ;
453     ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
454  
455     //Charged veto
456 //    Double_t cpvR=calo->GetEmcCpvDistance() ; 
457     Int_t itr=calo->GetNTracksMatched();  //number of track matched
458     ph->SetCPVBit(itr>0) ;  //Temporary cut !!!!
459
460     //Overall PID
461     Double_t pid[13];
462     calo->GetPID(pid);
463     ph->SetPCAPID(pid[AliAODCluster::kPhoton]>fPhotPID) ;
464       
465     //Distance to Bad
466     if(distBad>fMinDist3)
467       ph->SetDistToBad(2) ; 
468     else if(distBad>fMinDist2)
469         ph->SetDistToBad(1) ;
470     else
471        ph->SetDistToBad(0) ;
472  
473
474   }//loop
475   return kTRUE ;
476 }
477 //__________________________________________________________________
478 void  AliAnaPi0::FillHistograms() 
479 {
480   //Fill Real and Mixed invariant mass histograms
481   //then add current event to the buffer and 
482   //remove redundant events from buffer if necessary
483   
484   //Fill histograms only if list of particles for current event was created
485   if(fCurrentEvent==0)
486     return ;
487
488   //Fill Real distribution
489   Int_t nPhot = fCurrentEvent->GetEntriesFast() ;
490   for(Int_t i1=0; i1<nPhot-1; i1++){
491     AliCaloPhoton * p1 = static_cast<AliCaloPhoton*>(fCurrentEvent->At(i1)) ;
492     for(Int_t i2=i1+1; i2<nPhot; i2++){
493       AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(fCurrentEvent->At(i2)) ;
494       Double_t m =(*p1 + *p2).M() ;
495       Double_t pt=(*p1 + *p2).Pt();
496       Double_t a = TMath::Abs(p1->Energy()-p2->Energy())/(p1->Energy()+p2->Energy()) ;
497       for(Int_t ipid=0; ipid<fNPID; ipid++){
498         if(p1->IsPIDOK(ipid)&&p2->IsPIDOK(ipid)){   
499           fhRe1[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
500           if(p1->DistToBad()>0 && p2->DistToBad()>0){
501             fhRe2[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
502             if(p1->DistToBad()>1 && p2->DistToBad()>1){
503               fhRe3[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
504             }
505           }
506         }
507       } 
508     }
509   }
510
511   //Fill mixed
512   TList * evMixList=fEventsList[fCurCentrBin*fNZvertBin*fNrpBin+fCurZvertBin*fNrpBin+fCurRPBin] ;
513   Int_t nMixed = evMixList->GetSize() ;
514   for(Int_t ii=0; ii<nMixed; ii++){  
515     TClonesArray* ev2=dynamic_cast<TClonesArray*>(evMixList->At(ii));
516     Int_t nPhot2=ev2->GetEntriesFast() ;
517     for(Int_t i1=0; i1<nPhot; i1++){
518       AliCaloPhoton * p1 = static_cast<AliCaloPhoton*>(fCurrentEvent->At(i1)) ;
519       for(Int_t i2=0; i2<nPhot2; i2++){
520         AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(ev2->At(i2)) ;
521         Double_t m =(*p1 + *p2).M() ;
522         Double_t pt=(*p1 + *p2).Pt();
523         Double_t a = TMath::Abs(p1->Energy()-p2->Energy())/(p1->Energy()+p2->Energy()) ;
524         for(Int_t ipid=0; ipid<fNPID; ipid++){
525           if(p1->IsPIDOK(ipid)&&p2->IsPIDOK(ipid)){
526             fhMi1[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
527             if(p1->DistToBad()>0 && p2->DistToBad()>0){
528               fhMi2[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
529               if(p1->DistToBad()>1 && p2->DistToBad()>1){
530                 fhMi3[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
531               }
532             }
533           }
534         }
535       }
536     }
537   }
538
539   //Add current event to buffer and Remove redandant events 
540   if(fCurrentEvent->GetEntriesFast()>0){
541     evMixList->AddFirst(fCurrentEvent) ;
542     fCurrentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
543     if(evMixList->GetSize()>=fNmaxMixEv){
544       TClonesArray * tmp = dynamic_cast<TClonesArray*>(evMixList->Last()) ;
545       evMixList->RemoveLast() ;
546       delete tmp ;
547     }
548   } 
549   else{ //empty event
550    delete fCurrentEvent ;
551    fCurrentEvent=0 ; 
552   }
553 }
554 //______________________________________________________________________________
555 void AliAnaPi0::Terminate(Option_t *)
556 {
557 }