410d43803f5f047cf2efe9740ca3b55348503087
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / 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 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
23 //-- and Gustavo Conesa (INFN-Frascati)
24 //_________________________________________________________________________
25
26
27 // --- ROOT system ---
28 #include "TH3.h"
29 #include "TH2D.h"
30 //#include "Riostream.h"
31 #include "TCanvas.h"
32 #include "TPad.h"
33 #include "TROOT.h"
34 #include "TClonesArray.h"
35 #include "TObjString.h"
36
37 //---- AliRoot system ----
38 #include "AliAnaPi0.h"
39 #include "AliCaloTrackReader.h"
40 #include "AliCaloPID.h"
41 #include "AliStack.h"
42 #include "AliFiducialCut.h"
43 #include "TParticle.h"
44 #include "AliVEvent.h"
45 #include "AliESDCaloCluster.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliNeutralMesonSelection.h"
49 #include "AliMixedEvent.h"
50
51
52 ClassImp(AliAnaPi0)
53
54 //________________________________________________________________________________________________________________________________________________  
55 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
56 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
57 fNmaxMixEv(0), fCalorimeter(""),
58 fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), fMultiCutAna(kFALSE), 
59 fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fSameSM(kFALSE),
60 fhReMod(0x0),fhReDiffMod(0x0),
61 fhRe1(0x0),      fhMi1(0x0),      fhRe2(0x0),      fhMi2(0x0),      fhRe3(0x0),      fhMi3(0x0),
62 fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
63 fhRePtNCellAsymCuts(0x0), fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),  
64 fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
65 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
66 fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
67 {
68 //Default Ctor
69  InitParameters();
70  
71 }
72
73 //________________________________________________________________________________________________________________________________________________
74 AliAnaPi0::~AliAnaPi0() {
75   // Remove event containers
76   
77   if(fDoOwnMix && fEventsList){
78     for(Int_t ic=0; ic<fNCentrBin; ic++){
79       for(Int_t iz=0; iz<GetNZvertBin(); iz++){
80         for(Int_t irp=0; irp<GetNRPBin(); irp++){
81           fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
82           delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
83         }
84       }
85     }
86     delete[] fEventsList; 
87     fEventsList=0 ;
88   }
89         
90 }
91
92 //________________________________________________________________________________________________________________________________________________
93 void AliAnaPi0::InitParameters()
94 {
95 //Init parameters when first called the analysis
96 //Set default parameters
97   SetInputAODName("PWG4Particle");
98   
99   AddToHistogramsName("AnaPi0_");
100   fNModules = 12; // set maximum to maximum number of EMCAL modules
101   fNCentrBin = 1;
102 //  fNZvertBin = 1;
103 //  fNrpBin    = 1;
104   fNmaxMixEv = 10;
105  
106   fCalorimeter  = "PHOS";
107   fUseAngleCut = kFALSE;
108   
109   fMultiCutAna = kFALSE;
110   
111   fNPtCuts = 3;
112   fPtCuts[0] = 0.; fPtCuts[1] = 0.3;   fPtCuts[2] = 0.5;
113   for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
114   
115   fNAsymCuts = 4;
116   fAsymCuts[0] = 1.;  fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6;   fAsymCuts[3] = 0.1;    
117   for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
118
119   fNCellNCuts = 3;
120   fCellNCuts[0] = 0; fCellNCuts[1] = 1;   fCellNCuts[2] = 2;   
121   for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i]  = 0;
122
123   fNPIDBits = 2;
124   fPIDBits[0] = 0;   fPIDBits[1] = 2; //  fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut,  dispersion, neutral, dispersion&&neutral
125   for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
126
127 }
128
129
130 //________________________________________________________________________________________________________________________________________________
131 TObjString * AliAnaPi0::GetAnalysisCuts()
132 {  
133   //Save parameters used for analysis
134   TString parList ; //this will be list of parameters used for this analysis.
135   const Int_t buffersize = 255;
136   char onePar[buffersize] ;
137   snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
138   parList+=onePar ;     
139   snprintf(onePar,buffersize,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
140   parList+=onePar ;
141   snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
142   parList+=onePar ;
143   snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
144   parList+=onePar ;
145   snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
146   parList+=onePar ;
147   snprintf(onePar,buffersize,"Pair in same Module: %d \n",fSameSM) ;
148   parList+=onePar ;
149   snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
150   for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
151   parList+=onePar ;
152   snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
153   for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
154   parList+=onePar ;
155   snprintf(onePar,buffersize,"Cuts: \n") ;
156   parList+=onePar ;
157   snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
158   parList+=onePar ;
159   snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
160   parList+=onePar ;
161   snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
162   parList+=onePar ;
163   if(fMultiCutAna){
164     snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
165     for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
166     parList+=onePar ;
167     snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
168     for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
169     parList+=onePar ;
170   }
171   
172   return new TObjString(parList) ;      
173 }
174
175 //________________________________________________________________________________________________________________________________________________
176 TList * AliAnaPi0::GetCreateOutputObjects()
177 {  
178   // Create histograms to be saved in output file and 
179   // store them in fOutputContainer
180   
181   //create event containers
182   fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
183         
184   for(Int_t ic=0; ic<fNCentrBin; ic++){
185     for(Int_t iz=0; iz<GetNZvertBin(); iz++){
186       for(Int_t irp=0; irp<GetNRPBin(); irp++){
187         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
188         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
189       }
190     }
191   }
192   
193   TList * outputContainer = new TList() ; 
194   outputContainer->SetName(GetName()); 
195         
196   fhReMod     = new TH2D*[fNModules] ;
197   fhReDiffMod = new TH2D*[fNModules+1] ;
198   
199   fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
200   fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
201   fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
202   fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
203   fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
204   fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
205     
206   fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
207   fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
208   fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
209   fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
210   fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
211   fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
212     
213   const Int_t buffersize = 255;
214   char key[buffersize] ;
215   char title[buffersize] ;
216   
217   Int_t nptbins   = GetHistoPtBins();
218   Int_t nphibins  = GetHistoPhiBins();
219   Int_t netabins  = GetHistoEtaBins();
220   Float_t ptmax   = GetHistoPtMax();
221   Float_t phimax  = GetHistoPhiMax();
222   Float_t etamax  = GetHistoEtaMax();
223   Float_t ptmin   = GetHistoPtMin();
224   Float_t phimin  = GetHistoPhiMin();
225   Float_t etamin  = GetHistoEtaMin();   
226         
227   Int_t nmassbins = GetHistoMassBins();
228   Int_t nasymbins = GetHistoAsymmetryBins();
229   Float_t massmax = GetHistoMassMax();
230   Float_t asymmax = GetHistoAsymmetryMax();
231   Float_t massmin = GetHistoMassMin();
232   Float_t asymmin = GetHistoAsymmetryMin();
233   Int_t ntrmbins  = GetHistoTrackMultiplicityBins();
234   Int_t ntrmmax   = GetHistoTrackMultiplicityMax();
235   Int_t ntrmmin   = GetHistoTrackMultiplicityMin(); 
236
237   for(Int_t ic=0; ic<fNCentrBin; ic++){
238     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
239       for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
240         Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
241         //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
242         //Distance to bad module 1
243         snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
244         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
245                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
246         fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
247         fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
248         fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
249         //printf("name: %s\n ",fhRe1[index]->GetName());
250         outputContainer->Add(fhRe1[index]) ;
251         
252         //Distance to bad module 2
253         snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
254         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
255                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
256         fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
257         fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
258         fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
259         outputContainer->Add(fhRe2[index]) ;
260         
261         //Distance to bad module 3
262         snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
263         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
264                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
265         fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
266         fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
267         fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
268         outputContainer->Add(fhRe3[index]) ;
269         
270         //Inverse pT 
271         //Distance to bad module 1
272         snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
273         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
274                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
275         fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
276         fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
277         fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
278         outputContainer->Add(fhReInvPt1[index]) ;
279         
280         //Distance to bad module 2
281         snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
282         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
283                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
284         fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
285         fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
286         fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
287         outputContainer->Add(fhReInvPt2[index]) ;
288         
289         //Distance to bad module 3
290         snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
291         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
292                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
293         fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
294         fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
295         fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
296         outputContainer->Add(fhReInvPt3[index]) ;
297         
298         if(fDoOwnMix){
299           //Distance to bad module 1
300           snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
301           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
302                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
303           fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
304           fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
305           fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
306           outputContainer->Add(fhMi1[index]) ;
307           
308           //Distance to bad module 2
309           snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
310           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
311                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
312           fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
313           fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
314           fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
315           outputContainer->Add(fhMi2[index]) ;
316           
317           //Distance to bad module 3
318           snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
319           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
320                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
321           fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
322           fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
323           fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
324           outputContainer->Add(fhMi3[index]) ;
325           
326           //Inverse pT
327           //Distance to bad module 1
328           snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
329           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
330                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
331           fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
332           fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
333           fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
334           outputContainer->Add(fhMiInvPt1[index]) ;
335           
336           //Distance to bad module 2
337           snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
338           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
339                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
340           fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
341           fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
342           fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
343           outputContainer->Add(fhMiInvPt2[index]) ;
344           
345           //Distance to bad module 3
346           snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
347           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
348                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
349           fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
350           fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
351           fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
352           outputContainer->Add(fhMiInvPt3[index]) ;
353         } 
354       }
355     }
356   }
357   
358   fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
359   fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
360   fhRePtAsym->SetYTitle("Asymmetry");
361   outputContainer->Add(fhRePtAsym);
362   
363   fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
364   fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
365   fhRePtAsymPi0->SetYTitle("Asymmetry");
366   outputContainer->Add(fhRePtAsymPi0);
367
368   fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
369   fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
370   fhRePtAsymEta->SetYTitle("Asymmetry");
371   outputContainer->Add(fhRePtAsymEta);
372   
373   if(fMultiCutAna){
374     
375     fhRePIDBits         = new TH2D*[fNPIDBits];
376     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
377       snprintf(key,   buffersize,"hRe_pidbit%d",ipid) ;
378       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
379       fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
380       fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
381       fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
382       outputContainer->Add(fhRePIDBits[ipid]) ;
383     }// pid bit loop
384     
385     fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
386     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
387       for(Int_t icell=0; icell<fNCellNCuts; icell++){
388         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
389           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
390           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
391           Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
392           //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
393           fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
394           fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
395           fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
396           outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
397         }
398       }
399     }
400     
401     fhRePtMult = new TH3D*[fNAsymCuts] ;
402     for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
403       fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
404                                    nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
405       fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
406       fhRePtMult[iasym]->SetYTitle("Track multiplicity");
407       fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
408       outputContainer->Add(fhRePtMult[iasym]) ;
409     }
410     
411   }// multi cuts analysis
412   
413   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
414                     GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
415   outputContainer->Add(fhEvents) ;
416         
417   fhRealOpeningAngle  = new TH2D
418   ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5); 
419   fhRealOpeningAngle->SetYTitle("#theta(rad)");
420   fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
421   outputContainer->Add(fhRealOpeningAngle) ;
422   
423   fhRealCosOpeningAngle  = new TH2D
424   ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1); 
425   fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
426   fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
427   outputContainer->Add(fhRealCosOpeningAngle) ;
428         
429   //Histograms filled only if MC data is requested      
430   if(IsDataMC()){
431
432     fhPrimPt     = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
433     fhPrimAccPt  = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
434     outputContainer->Add(fhPrimPt) ;
435     outputContainer->Add(fhPrimAccPt) ;
436     
437     fhPrimY      = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
438     outputContainer->Add(fhPrimY) ;
439     
440     fhPrimAccY   = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
441     outputContainer->Add(fhPrimAccY) ;
442     
443     fhPrimPhi    = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
444     outputContainer->Add(fhPrimPhi) ;
445     
446     fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
447     outputContainer->Add(fhPrimAccPhi) ;
448     
449     
450     fhPrimOpeningAngle  = new TH2D
451     ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5); 
452     fhPrimOpeningAngle->SetYTitle("#theta(rad)");
453     fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
454     outputContainer->Add(fhPrimOpeningAngle) ;
455     
456     fhPrimCosOpeningAngle  = new TH2D
457     ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1); 
458     fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
459     fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
460     outputContainer->Add(fhPrimCosOpeningAngle) ;
461     
462   }
463   
464   TString * pairname = new TString[fNModules];
465   if(fCalorimeter=="EMCAL"){
466     pairname[0]="A side (0-2)"; 
467     pairname[1]="C side (1-3)";
468     pairname[2]="Sector 0 (0-1)"; 
469     pairname[3]="Sector 1 (2-3)";
470     for(Int_t i = 4 ; i < fNModules ; i++) pairname[i]="";}
471   if(fCalorimeter=="PHOS") {
472     pairname[0]="(0-1)"; 
473     pairname[1]="(0-2)";
474     pairname[2]="(1-2)";
475     for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
476
477   for(Int_t imod=0; imod<fNModules; imod++){
478     //Module dependent invariant mass
479     snprintf(key, buffersize,"hReMod_%d",imod) ;
480     snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
481     fhReMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
482     fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
483     fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
484     outputContainer->Add(fhReMod[imod]) ;
485
486     snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
487     snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
488     fhReDiffMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
489     fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
490     fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
491     outputContainer->Add(fhReDiffMod[imod]) ;
492   }
493   
494   delete [] pairname;
495   
496   snprintf(key, buffersize,"hReDiffMod_%d",fNModules) ;
497   snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for all Modules Combination") ;
498   fhReDiffMod[fNModules]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
499   outputContainer->Add(fhReDiffMod[fNModules]) ;
500   
501   
502 //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
503 //  
504 //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
505 //  
506 //  }
507   
508   return outputContainer;
509 }
510
511 //_________________________________________________________________________________________________________________________________________________
512 void AliAnaPi0::Print(const Option_t * /*opt*/) const
513 {
514   //Print some relevant parameters set for the analysis
515   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
516   AliAnaPartCorrBaseClass::Print(" ");
517
518   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
519   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
520   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
521   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
522   printf("Pair in same Module: %d \n",fSameSM) ;
523   printf("Cuts: \n") ;
524   printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
525   printf("Number of modules:             %d \n",fNModules) ;
526   printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
527   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
528   printf("\tasymmetry < ");
529   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
530   printf("\n");
531   
532   printf("PID selection bits: n = %d, \n",fNPIDBits) ;
533   printf("\tPID bit = ");
534   for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
535   printf("\n");
536   
537   if(fMultiCutAna){
538     printf("pT cuts: n = %d, \n",fNPtCuts) ;
539     printf("\tpT > ");
540     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
541     printf("GeV/c\n");
542     
543     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
544     printf("\tnCell > ");
545     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
546     printf("\n");
547
548   }
549   printf("------------------------------------------------------\n") ;
550
551
552 //_____________________________________________________________
553 void AliAnaPi0::FillAcceptanceHistograms(){
554   //Fill acceptance histograms if MC data is available
555   
556   if(IsDataMC() && GetReader()->ReadStack()){   
557     AliStack * stack = GetMCStack();
558     if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
559       for(Int_t i=0 ; i<stack->GetNprimary(); i++){
560         TParticle * prim = stack->Particle(i) ;
561         if(prim->GetPdgCode() == 111){
562           Double_t pi0Pt = prim->Pt() ;
563           //printf("pi0, pt %2.2f\n",pi0Pt);
564           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception          
565           Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
566           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
567           if(TMath::Abs(pi0Y) < 0.5){
568             fhPrimPt->Fill(pi0Pt) ;
569           }
570           fhPrimY  ->Fill(pi0Y) ;
571           fhPrimPhi->Fill(phi) ;
572           
573           //Check if both photons hit Calorimeter
574           Int_t iphot1=prim->GetFirstDaughter() ;
575           Int_t iphot2=prim->GetLastDaughter() ;
576           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
577             TParticle * phot1 = stack->Particle(iphot1) ;
578             TParticle * phot2 = stack->Particle(iphot2) ;
579             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
580               //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
581               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
582               
583               TLorentzVector lv1, lv2;
584               phot1->Momentum(lv1);
585               phot2->Momentum(lv2);
586               
587               Bool_t inacceptance = kFALSE;
588               if(fCalorimeter == "PHOS"){
589                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
590                   Int_t mod ;
591                   Double_t x,z ;
592                   if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x)) 
593                     inacceptance = kTRUE;
594                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
595                 }
596                 else{
597                   
598                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter)) 
599                     inacceptance = kTRUE ;
600                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
601                 }
602                 
603               }    
604               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
605                 if(GetEMCALGeometry()){
606                   if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2)) 
607                     inacceptance = kTRUE;
608                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
609                 }
610                 else{
611                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter)) 
612                     inacceptance = kTRUE ;
613                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
614                 }
615               }   
616               
617               if(inacceptance){
618                 
619                 fhPrimAccPt->Fill(pi0Pt) ;
620                 fhPrimAccPhi->Fill(phi) ;
621                 fhPrimAccY->Fill(pi0Y) ;
622                 Double_t angle  = lv1.Angle(lv2.Vect());
623                 fhPrimOpeningAngle   ->Fill(pi0Pt,angle);
624                 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
625                 
626               }//Accepted
627             }// 2 photons      
628           }//Check daughters exist
629         }// Primary pi0
630       }//loop on primaries      
631     }//stack exists and data is MC
632   }//read stack
633   else if(GetReader()->ReadAODMCParticles()){
634     if(GetDebug() >= 0)  printf("AliAnaPi0::FillAcceptanceHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
635   }     
636 }
637
638 //____________________________________________________________________________________________________________________________________________________
639 void AliAnaPi0::MakeAnalysisFillHistograms() 
640 {
641   //Process one event and extract photons from AOD branch 
642   // filled with AliAnaPhoton and fill histos with invariant mass
643   
644   //In case of MC data, fill acceptance histograms
645   FillAcceptanceHistograms();
646   
647   //Apply some cuts on event: vertex position and centrality range  
648   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
649   if(IsBadRun(iRun)) return ;   
650   
651   Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
652   if(GetDebug() > 1) 
653     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
654   if(nPhot < 2 )
655     return ; 
656   Int_t module1 = -1;
657   Int_t module2 = -1;
658   Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex 
659   Int_t evtIndex1 = 0 ; 
660   Int_t currentEvtIndex = -1 ; 
661   Int_t curCentrBin = 0 ; 
662   Int_t curRPBin    = 0 ; 
663   Int_t curZvertBin = 0 ;
664   
665   for(Int_t i1=0; i1<nPhot-1; i1++){
666     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
667     // get the event index in the mixed buffer where the photon comes from 
668     // in case of mixing with analysis frame, not own mixing
669     evtIndex1 = GetEventIndex(p1, vert) ; 
670     //printf("charge = %d\n", track->Charge());
671     if ( evtIndex1 == -1 )
672       return ; 
673     if ( evtIndex1 == -2 )
674       continue ; 
675     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
676     if (evtIndex1 != currentEvtIndex) {
677       curCentrBin = GetEventCentrality();
678       curRPBin    = 0 ;
679       curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
680       fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
681       currentEvtIndex = evtIndex1 ; 
682       //if(GetDebug() > 1) 
683         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d\n",curCentrBin,curRPBin,curZvertBin);
684     }
685     
686     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
687     
688     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
689     //Get Module number
690     module1 = GetModuleNumber(p1);
691     for(Int_t i2=i1+1; i2<nPhot; i2++){
692       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
693       Int_t evtIndex2 = GetEventIndex(p2, vert) ; 
694       if ( evtIndex2 == -1 )
695         return ; 
696       if ( evtIndex2 == -2 )
697         continue ;    
698       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
699         continue ;
700       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
701       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
702       //Get module number
703       module2 = GetModuleNumber(p2);
704       Double_t m  = (photon1 + photon2).M() ;
705       Double_t pt = (photon1 + photon2).Pt();
706       Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
707       if(GetDebug() > 2)
708         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
709                p1->Pt(), p2->Pt(), pt,m,a);
710       //Check if opening angle is too large or too small compared to what is expected   
711       Double_t angle   = photon1.Angle(photon2.Vect());
712       //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
713       //printf("angle %f\n",angle);
714       if(fUseAngleCut && angle < 0.1) 
715         continue;
716
717       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
718       if(a < fAsymCuts[0]){
719         if(module1==module2 && module1 >=0 && module1<fNModules)
720           fhReMod[module1]->Fill(pt,m) ;
721         else  
722           fhReDiffMod[fNModules]->Fill(pt,m) ;
723         
724         if(fCalorimeter=="EMCAL"){
725           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ; 
726           if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,m) ; 
727           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,m) ;
728           if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,m) ; 
729         }
730         else {
731           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ; 
732           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[1]->Fill(pt,m) ; 
733           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffMod[2]->Fill(pt,m) ;
734         }
735       }
736       
737       //In case we want only pairs in same (super) module, check their origin.
738       Bool_t ok = kTRUE;
739       if(fSameSM && module1!=module2) ok=kFALSE;
740       if(ok){
741         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
742         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
743           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){ 
744             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
745               if(a < fAsymCuts[iasym]){
746                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
747                 //printf("cen %d, pid %d, asy %d, Index %d\n",curCentrBin,ipid,iasym,index);
748                 fhRe1     [index]->Fill(pt,m);
749                 fhReInvPt1[index]->Fill(pt,m,1./pt) ;
750                 if(p1->DistToBad()>0 && p2->DistToBad()>0){
751                   fhRe2     [index]->Fill(pt,m) ;
752                   fhReInvPt2[index]->Fill(pt,m,1./pt) ;
753                   if(p1->DistToBad()>1 && p2->DistToBad()>1){
754                     fhRe3     [index]->Fill(pt,m) ;
755                     fhReInvPt3[index]->Fill(pt,m,1./pt) ;
756                   }//assymetry cut
757                 }// asymmetry cut loop
758               }// bad 3
759             }// bad2
760           }// bad 1
761         }// pid bit loop
762         
763         //Fill histograms with opening angle
764         fhRealOpeningAngle   ->Fill(pt,angle);
765         fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
766         
767         //Fill histograms with pair assymmetry
768         fhRePtAsym->Fill(pt,a);
769         if(m > 0.10 && m < 0.16) fhRePtAsymPi0->Fill(pt,a);
770         if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
771         
772         //Multi cuts analysis 
773         if(fMultiCutAna){
774           //Histograms for different PID bits selection
775           for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
776             
777             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    && 
778                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
779             
780             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
781           } // pid bit cut loop
782           
783           //Several pt,ncell and asymmetry cuts
784           //Get the number of cells
785           Int_t ncell1 = 0;
786           Int_t ncell2 = 0;
787           AliVEvent * event = GetReader()->GetInputEvent();
788           if(event){
789             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
790               AliVCluster *cluster = event->GetCaloCluster(iclus);
791               
792               Bool_t is = kFALSE;
793               if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
794               else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
795               
796               if(is){
797                 if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
798                 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
799               } // PHOS or EMCAL cluster as requested in analysis
800               
801               if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
802               
803             }
804             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
805           }
806           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){          
807             for(Int_t icell=0; icell<fNCellNCuts; icell++){
808               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
809                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
810                 if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        && 
811                    a        <   fAsymCuts[iasym]                                    && 
812                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]) fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
813                 
814                 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
815               }// pid bit cut loop
816             }// icell loop
817           }// pt cut loop
818           for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
819             if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
820           }
821           
822         }// multiple cuts analysis
823       }// ok if same sm
824     }// second same event particle
825   }// first cluster
826   
827   if(fDoOwnMix){
828     //Fill mixed
829     TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
830     Int_t nMixed = evMixList->GetSize() ;
831     for(Int_t ii=0; ii<nMixed; ii++){  
832       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
833       Int_t nPhot2=ev2->GetEntriesFast() ;
834       Double_t m = -999;
835       if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
836       
837       for(Int_t i1=0; i1<nPhot; i1++){
838         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
839         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
840         module1 = GetModuleNumber(p1);
841         for(Int_t i2=0; i2<nPhot2; i2++){
842           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
843           
844           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
845           m =           (photon1+photon2).M() ; 
846           Double_t pt = (photon1 + photon2).Pt();
847           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
848           
849           //Check if opening angle is too large or too small compared to what is expected
850           Double_t angle   = photon1.Angle(photon2.Vect());
851           //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
852           if(fUseAngleCut && angle < 0.1) continue;  
853           
854           if(GetDebug() > 2)
855             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
856                    p1->Pt(), p2->Pt(), pt,m,a); 
857           //In case we want only pairs in same (super) module, check their origin.
858           module2 = GetModuleNumber(p2);
859           Bool_t ok = kTRUE;
860           if(fSameSM && module1!=module2) ok=kFALSE;
861           if(ok){
862             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){ 
863               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
864                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
865                   if(a < fAsymCuts[iasym]){
866                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
867                     fhMi1     [index]->Fill(pt,m) ;
868                     fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
869                     if(p1->DistToBad()>0 && p2->DistToBad()>0){
870                       fhMi2     [index]->Fill(pt,m) ;
871                       fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
872                       if(p1->DistToBad()>1 && p2->DistToBad()>1){
873                         fhMi3     [index]->Fill(pt,m) ;
874                         fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
875                       }
876                     }
877                   }//Asymmetry cut
878                 }// Asymmetry loop
879               }//PID cut
880             }//loop for histograms
881           }//ok
882         }// second cluster loop
883       }//first cluster loop
884     }//loop on mixed events
885     
886     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
887     //Add current event to buffer and Remove redundant events 
888     if(currentEvent->GetEntriesFast()>0){
889       evMixList->AddFirst(currentEvent) ;
890       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
891       if(evMixList->GetSize()>=fNmaxMixEv)
892       {
893         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
894         evMixList->RemoveLast() ;
895         delete tmp ;
896       }
897     } 
898     else{ //empty event
899       delete currentEvent ;
900       currentEvent=0 ; 
901     }
902   }// DoOwnMix
903   
904 }       
905
906 //________________________________________________________________________
907 void AliAnaPi0::ReadHistograms(TList* outputList)
908 {
909   // Needed when Terminate is executed in distributed environment
910   // Refill analysis histograms of this class with corresponding histograms in output list. 
911   
912   // Histograms of this analsys are kept in the same list as other analysis, recover the position of
913   // the first one and then add the next.
914   Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
915   
916   if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
917   if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
918   if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
919   if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
920   if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
921   if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;       
922   if(!fhReInvPt1) fhReInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
923   if(!fhReInvPt2) fhReInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
924   if(!fhReInvPt3) fhReInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
925   if(!fhMiInvPt1) fhMiInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
926   if(!fhMiInvPt2) fhMiInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
927   if(!fhMiInvPt3) fhMiInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;    
928   if(!fhReMod)    fhReMod     = new TH2D*[fNModules]   ;        
929   if(!fhReDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;        
930
931   for(Int_t ic=0; ic<fNCentrBin; ic++){
932     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
933       for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
934         Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
935
936         fhRe1[ihisto] = (TH2D*) outputList->At(index++);
937         fhRe2[ihisto] = (TH2D*) outputList->At(index++);
938         fhRe3[ihisto] = (TH2D*) outputList->At(index++);
939       
940         fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
941         fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
942         fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
943       
944         if(fDoOwnMix){
945           fhMi1[ihisto] = (TH2D*) outputList->At(index++);
946           fhMi2[ihisto] = (TH2D*) outputList->At(index++);
947           fhMi3[ihisto] = (TH2D*) outputList->At(index++);
948       
949           fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
950           fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
951           fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++); 
952         }//Own mix
953       }//asymmetry loop
954     }// pid loop
955   }// centrality loop
956   
957   fhRePtAsym    = (TH2D*)outputList->At(index++);
958   fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
959   fhRePtAsymEta = (TH2D*)outputList->At(index++);
960   
961   if(fMultiCutAna){
962     
963     if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
964     if(!fhRePIDBits)         fhRePIDBits         = new TH2D*[fNPIDBits];
965
966     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
967       fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
968     }// ipid loop
969     
970     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
971       for(Int_t icell=0; icell<fNCellNCuts; icell++){
972         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
973           fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
974         }// iasym
975       }// icell loop
976     }// ipt loop
977     
978     if(!fhRePtMult) fhRePtMult  = new TH3D*[fNAsymCuts]  ;
979     for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
980       fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
981   }// multi cut analysis 
982   
983   fhEvents = (TH3D *) outputList->At(index++); 
984   
985   fhRealOpeningAngle     = (TH2D*)  outputList->At(index++);
986   fhRealCosOpeningAngle  = (TH2D*)  outputList->At(index++);
987   
988   //Histograms filled only if MC data is requested      
989   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
990     fhPrimPt     = (TH1D*)  outputList->At(index++);
991     fhPrimAccPt  = (TH1D*)  outputList->At(index++);
992     fhPrimY      = (TH1D*)  outputList->At(index++);
993     fhPrimAccY   = (TH1D*)  outputList->At(index++);
994     fhPrimPhi    = (TH1D*)  outputList->At(index++);
995     fhPrimAccPhi = (TH1D*)  outputList->At(index++);
996   }
997   
998   for(Int_t imod=0; imod < fNModules; imod++)
999     fhReMod[imod] = (TH2D*) outputList->At(index++);
1000   
1001   
1002 }
1003
1004
1005 //____________________________________________________________________________________________________________________________________________________
1006 void AliAnaPi0::Terminate(TList* outputList) 
1007 {
1008   //Do some calculations and plots from the final histograms.
1009   
1010   printf(" *** %s Terminate:\n", GetName()) ; 
1011   
1012   //Recover histograms from output histograms list, needed for distributed analysis.    
1013   ReadHistograms(outputList);
1014   
1015   if (!fhRe1) {
1016     printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
1017     return;
1018   }
1019   
1020   printf("AliAnaPi0::Terminate()         Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
1021     
1022   const Int_t buffersize = 255;
1023
1024   char nameIM[buffersize];
1025   snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
1026   TCanvas  * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
1027   cIM->Divide(2, 2);
1028   
1029   cIM->cd(1) ; 
1030   //gPad->SetLogy();
1031   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
1032   hIMAllPt->SetLineColor(2);
1033   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
1034   hIMAllPt->Draw();
1035
1036   cIM->cd(2) ; 
1037   TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
1038 //  hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
1039 //  TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
1040   hIMPt5->SetLineColor(2);  
1041   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
1042   hIMPt5->Draw();
1043   
1044   cIM->cd(3) ; 
1045   TH1D * hIMPt10 =  (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
1046 //  hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
1047 //  TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
1048   hIMPt10->SetLineColor(2);  
1049   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
1050   hIMPt10->Draw();
1051   
1052   cIM->cd(4) ; 
1053   TH1D * hIMPt20 =  (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
1054  // TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
1055 //  hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
1056 //  TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
1057   hIMPt20->SetLineColor(2);  
1058   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
1059   hIMPt20->Draw();
1060    
1061   char nameIMF[buffersize];
1062   snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
1063   cIM->Print(nameIMF);
1064
1065   char namePt[buffersize];
1066   snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
1067   TCanvas  * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
1068   cPt->Divide(2, 2);
1069
1070   cPt->cd(1) ; 
1071   //gPad->SetLogy();
1072   TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
1073   hPt->SetLineColor(2);
1074   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
1075   hPt->Draw();
1076
1077   cPt->cd(2) ; 
1078   TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21)); 
1079 //  TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
1080 //  hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
1081 //  TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
1082   hPtIM1->SetLineColor(2);  
1083   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
1084   hPtIM1->Draw();
1085   
1086   cPt->cd(3) ; 
1087   TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17)); 
1088 //  TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
1089 //  hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
1090 //  TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
1091   hPtIM2->SetLineColor(2);  
1092   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
1093   hPtIM2->Draw();
1094
1095   cPt->cd(4) ; 
1096   TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15)); 
1097 //  TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
1098 //  hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
1099 //  TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
1100   hPtIM3->SetLineColor(2);  
1101   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
1102   hPtIM3->Draw();
1103    
1104   char namePtF[buffersize];
1105   snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
1106   cPt->Print(namePtF);
1107
1108   char line[buffersize] ; 
1109   snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ; 
1110   gROOT->ProcessLine(line);
1111   snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data()); 
1112   gROOT->ProcessLine(line);
1113  
1114   printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1115
1116 }
1117   //____________________________________________________________________________________________________________________________________________________
1118 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)  
1119 {
1120   // retieves the event index and checks the vertex
1121   //    in the mixed buffer returns -2 if vertex NOK
1122   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
1123   
1124   Int_t evtIndex = -1 ; 
1125   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
1126     
1127     if (GetMixedEvent()){
1128       
1129       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
1130       GetVertex(vert,evtIndex); 
1131       
1132       if(TMath::Abs(vert[2])> GetZvertexCut())
1133         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1134     } else {// Single event
1135       
1136       GetVertex(vert);
1137       
1138       if(TMath::Abs(vert[2])> GetZvertexCut())
1139         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1140       else 
1141         evtIndex = 0 ;
1142     }
1143   }//No MC reader
1144   else {
1145     evtIndex = 0;
1146     vert[0] = 0. ; 
1147     vert[1] = 0. ; 
1148     vert[2] = 0. ; 
1149   }
1150   
1151   return evtIndex ; 
1152 }
1153