change all printf's to AliDebug/AliInfo/AliWarning
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / 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
16 //_________________________________________________________________________
17 // Class to collect two-photon invariant mass distributions for
18 // extracting raw pi0 yield.
19 // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
20 // it will do nothing if executed alone
21 //
22 //-- Author: Dmitri Peressounko (RRC "KI")
23 //-- Adapted to CaloTrackCorr frame by Lamia Benhabib (SUBATECH)
24 //-- and Gustavo Conesa (INFN-Frascati)
25 //_________________________________________________________________________
26
27
28 // --- ROOT system ---
29 #include "TH3.h"
30 #include "TH2F.h"
31 //#include "Riostream.h"
32 #include "TCanvas.h"
33 #include "TPad.h"
34 #include "TROOT.h"
35 #include "TClonesArray.h"
36 #include "TObjString.h"
37 #include "TDatabasePDG.h"
38
39 //---- AliRoot system ----
40 #include "AliAnaPi0.h"
41 #include "AliCaloTrackReader.h"
42 #include "AliCaloPID.h"
43 #include "AliStack.h"
44 #include "AliFiducialCut.h"
45 #include "TParticle.h"
46 #include "AliVEvent.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliNeutralMesonSelection.h"
51 #include "AliMixedEvent.h"
52 #include "AliAODMCParticle.h"
53
54 // --- Detectors ---
55 #include "AliPHOSGeoUtils.h"
56 #include "AliEMCALGeometry.h"
57
58 ClassImp(AliAnaPi0)
59
60 //______________________________________________________
61 AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
62 fEventsList(0x0),
63 fNModules(22),
64 fUseAngleCut(kFALSE),        fUseAngleEDepCut(kFALSE),     fAngleCut(0),                 fAngleMaxCut(7.),
65 fMultiCutAna(kFALSE),        fMultiCutAnaSim(kFALSE),
66 fNPtCuts(0),                 fNAsymCuts(0),                fNCellNCuts(0),               fNPIDBits(0),
67 fMakeInvPtPlots(kFALSE),     fSameSM(kFALSE),
68 fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
69 fFillBadDistHisto(kFALSE),   fFillSSCombinations(kFALSE),
70 fFillAngleHisto(kFALSE),     fFillAsymmetryHisto(kFALSE),  fFillOriginHisto(0),          fFillArmenterosThetaStar(0),
71 fCheckAccInSector(kFALSE),
72 fPhotonMom1(),               fPhotonMom1Boost(),           fPhotonMom2(),                fPi0Mom(),
73 fProdVertex(),
74
75 //Histograms
76 fhAverTotECluster(0),        fhAverTotECell(0),            fhAverTotECellvsCluster(0),
77 fhEDensityCluster(0),        fhEDensityCell(0),            fhEDensityCellvsCluster(0),
78 fhReMod(0x0),                fhReSameSideEMCALMod(0x0),    fhReSameSectorEMCALMod(0x0),  fhReDiffPHOSMod(0x0),
79 fhMiMod(0x0),                fhMiSameSideEMCALMod(0x0),    fhMiSameSectorEMCALMod(0x0),  fhMiDiffPHOSMod(0x0),
80 fhReConv(0x0),               fhMiConv(0x0),                fhReConv2(0x0),  fhMiConv2(0x0),
81 fhRe1(0x0),                  fhMi1(0x0),                   fhRe2(0x0),                   fhMi2(0x0),
82 fhRe3(0x0),                  fhMi3(0x0),                   fhReInvPt1(0x0),              fhMiInvPt1(0x0),
83 fhReInvPt2(0x0),             fhMiInvPt2(0x0),              fhReInvPt3(0x0),              fhMiInvPt3(0x0),
84 fhRePtNCellAsymCuts(0x0),    fhMiPtNCellAsymCuts(0x0),     fhRePtNCellAsymCutsSM(),
85 fhRePIDBits(0x0),            fhRePtMult(0x0),              fhReSS(),
86 fhRePtAsym(0x0),             fhRePtAsymPi0(0x0),           fhRePtAsymEta(0x0),
87 fhEventBin(0),               fhEventMixBin(0),
88 fhCentrality(0x0),           fhCentralityNoPair(0x0),
89 fhEventPlaneResolution(0x0),
90 fhRealOpeningAngle(0x0),     fhRealCosOpeningAngle(0x0),   fhMixedOpeningAngle(0x0),     fhMixedCosOpeningAngle(0x0),
91 // MC histograms
92 fhPrimPi0E(0x0),             fhPrimPi0Pt(0x0),
93 fhPrimPi0AccE(0x0),          fhPrimPi0AccPt(0x0),
94 fhPrimPi0Y(0x0),             fhPrimPi0AccY(0x0),
95 fhPrimPi0Yeta(0x0),          fhPrimPi0YetaYcut(0x0),       fhPrimPi0AccYeta(0x0),
96 fhPrimPi0Phi(0x0),           fhPrimPi0AccPhi(0x0),
97 fhPrimPi0OpeningAngle(0x0),  fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
98 fhPrimPi0PtCentrality(0),    fhPrimPi0PtEventPlane(0),
99 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
100 fhPrimEtaE(0x0),             fhPrimEtaPt(0x0),
101 fhPrimEtaAccE(0x0),          fhPrimEtaAccPt(0x0),
102 fhPrimEtaY(0x0),             fhPrimEtaAccY(0x0),
103 fhPrimEtaYeta(0x0),          fhPrimEtaYetaYcut(0x0),       fhPrimEtaAccYeta(0x0),
104 fhPrimEtaPhi(0x0),           fhPrimEtaAccPhi(0x0),
105 fhPrimEtaOpeningAngle(0x0),  fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
106 fhPrimEtaPtCentrality(0),    fhPrimEtaPtEventPlane(0),
107 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
108 fhPrimPi0PtOrigin(0x0),      fhPrimEtaPtOrigin(0x0),
109 fhMCOrgMass(),               fhMCOrgAsym(),                fhMCOrgDeltaEta(),            fhMCOrgDeltaPhi(),
110 fhMCPi0MassPtRec(),          fhMCPi0MassPtTrue(),          fhMCPi0PtTruePtRec(),
111 fhMCEtaMassPtRec(),          fhMCEtaMassPtTrue(),          fhMCEtaPtTruePtRec(),
112 fhMCPi0PtOrigin(0x0),        fhMCEtaPtOrigin(0x0),
113 fhMCPi0ProdVertex(0),        fhMCEtaProdVertex(0),
114 fhPrimPi0ProdVertex(0),      fhPrimEtaProdVertex(0),
115 fhReMCFromConversion(0),     fhReMCFromNotConversion(0),   fhReMCFromMixConversion(0),
116 fhCosThStarPrimPi0(0),       fhCosThStarPrimEta(0)//,
117 {
118   //Default Ctor
119   
120   InitParameters();
121   
122   for(Int_t i = 0; i < 4; i++)
123   {
124     fhArmPrimEta[i] = 0;
125     fhArmPrimPi0[i] = 0;
126   }
127 }
128
129 //_____________________
130 AliAnaPi0::~AliAnaPi0()
131 {
132   // Remove event containers
133   
134   if(DoOwnMix() && fEventsList)
135   {
136     for(Int_t ic=0; ic<GetNCentrBin(); ic++)
137     {
138       for(Int_t iz=0; iz<GetNZvertBin(); iz++)
139       {
140         for(Int_t irp=0; irp<GetNRPBin(); irp++)
141         {
142           Int_t bin = GetEventMixBin(ic,iz,irp);
143           fEventsList[bin]->Delete() ;
144           delete fEventsList[bin] ;
145         }
146       }
147     }
148     delete[] fEventsList;
149   }
150   
151 }
152
153 //______________________________
154 void AliAnaPi0::InitParameters()
155 {
156   //Init parameters when first called the analysis
157   //Set default parameters
158   SetInputAODName("PWG4Particle");
159   
160   AddToHistogramsName("AnaPi0_");
161   
162   fUseAngleCut = kFALSE;
163   fUseAngleEDepCut = kFALSE;
164   fAngleCut    = 0.;
165   fAngleMaxCut = TMath::Pi();
166   
167   fMultiCutAna = kFALSE;
168   
169   fNPtCuts = 3;
170   fPtCuts[0] = 0.; fPtCuts[1] = 0.3;   fPtCuts[2] = 0.5;
171   for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
172   
173   fNAsymCuts = 2;
174   fAsymCuts[0] = 1.;  fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; //  fAsymCuts[3] = 0.1;
175   for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
176   
177   fNCellNCuts = 3;
178   fCellNCuts[0] = 0; fCellNCuts[1] = 1;   fCellNCuts[2] = 2;
179   for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i]  = 0;
180   
181   fNPIDBits = 2;
182   fPIDBits[0] = 0;   fPIDBits[1] = 2; //  fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut,  dispersion, neutral, dispersion&&neutral
183   for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
184   
185 }
186
187
188 //_______________________________________
189 TObjString * AliAnaPi0::GetAnalysisCuts()
190 {
191   //Save parameters used for analysis
192   TString parList ; //this will be list of parameters used for this analysis.
193   const Int_t buffersize = 255;
194   char onePar[buffersize] ;
195   snprintf(onePar,buffersize,"--- AliAnaPi0 ---:") ;
196   parList+=onePar ;
197   snprintf(onePar,buffersize,"Number of bins in Centrality:  %d;",GetNCentrBin()) ;
198   parList+=onePar ;
199   snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d;",GetNZvertBin()) ;
200   parList+=onePar ;
201   snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d;",GetNRPBin()) ;
202   parList+=onePar ;
203   snprintf(onePar,buffersize,"Depth of event buffer: %d;",GetNMaxEvMix()) ;
204   parList+=onePar ;
205   snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f;",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
206   parList+=onePar ;
207   snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
208   for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
209   parList+=onePar ;
210   snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =",fNPIDBits) ;
211   for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
212   parList+=onePar ;
213   snprintf(onePar,buffersize,"Cuts:") ;
214   parList+=onePar ;
215   snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f;",GetZvertexCut(),GetZvertexCut()) ;
216   parList+=onePar ;
217   snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
218   parList+=onePar ;
219   snprintf(onePar,buffersize,"Number of modules: %d:",fNModules) ;
220   parList+=onePar ;
221   if(fMultiCutAna)
222   {
223     snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
224     for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
225     parList+=onePar ;
226     snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
227     for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
228     parList+=onePar ;
229   }
230   
231   return new TObjString(parList) ;
232 }
233
234 //_________________________________________
235 TList * AliAnaPi0::GetCreateOutputObjects()
236 {
237   // Create histograms to be saved in output file and
238   // store them in fOutputContainer
239   
240   // Init the number of modules, set in the class AliCalorimeterUtils
241   fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
242   if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
243   
244   //create event containers
245   fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
246   
247   for(Int_t ic=0; ic<GetNCentrBin(); ic++)
248   {
249     for(Int_t iz=0; iz<GetNZvertBin(); iz++)
250     {
251       for(Int_t irp=0; irp<GetNRPBin(); irp++)
252       {
253         Int_t bin = GetEventMixBin(ic,iz,irp);
254         fEventsList[bin] = new TList() ;
255         fEventsList[bin]->SetOwner(kFALSE);
256       }
257     }
258   }
259   
260   TList * outputContainer = new TList() ;
261   outputContainer->SetName(GetName());
262   
263   fhReMod                = new TH2F*[fNModules]   ;
264   fhMiMod                = new TH2F*[fNModules]   ;
265   
266   if(GetCalorimeter() == kPHOS)
267   {
268     fhReDiffPHOSMod        = new TH2F*[fNModules]   ;
269     fhMiDiffPHOSMod        = new TH2F*[fNModules]   ;
270   }
271   else
272   {
273     fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
274     fhReSameSideEMCALMod   = new TH2F*[fNModules-2] ;
275     fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
276     fhMiSameSideEMCALMod   = new TH2F*[fNModules-2] ;
277   }
278   
279   
280   fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
281   fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
282   if(fFillBadDistHisto)
283   {
284     fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
285     fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
286     fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
287     fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
288   }
289   if(fMakeInvPtPlots)
290   {
291     fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
292     fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
293     if(fFillBadDistHisto){
294       fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
295       fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
296       fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
297       fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
298     }
299   }
300   
301   const Int_t buffersize = 255;
302   char key[buffersize] ;
303   char title[buffersize] ;
304   
305   Int_t nptbins   = GetHistogramRanges()->GetHistoPtBins();
306   Int_t nphibins  = GetHistogramRanges()->GetHistoPhiBins();
307   Int_t netabins  = GetHistogramRanges()->GetHistoEtaBins();
308   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();
309   Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();
310   Float_t etamax  = GetHistogramRanges()->GetHistoEtaMax();
311   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();
312   Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();
313   Float_t etamin  = GetHistogramRanges()->GetHistoEtaMin();
314   
315   Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
316   Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins();
317   Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
318   Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax();
319   Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
320   Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
321   Int_t ntrmbins  = GetHistogramRanges()->GetHistoTrackMultiplicityBins();
322   Int_t ntrmmax   = GetHistogramRanges()->GetHistoTrackMultiplicityMax();
323   Int_t ntrmmin   = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
324   
325   if(fCheckConversion)
326   {
327     fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
328     fhReConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
329     fhReConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
330     outputContainer->Add(fhReConv) ;
331     
332     fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
333     fhReConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
334     fhReConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
335     outputContainer->Add(fhReConv2) ;
336     
337     if(DoOwnMix())
338     {
339       fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
340       fhMiConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
341       fhMiConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
342       outputContainer->Add(fhMiConv) ;
343       
344       fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
345       fhMiConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
346       fhMiConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
347       outputContainer->Add(fhMiConv2) ;
348     }
349   }
350   
351   for(Int_t ic=0; ic<GetNCentrBin(); ic++)
352   {
353     for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
354     {
355       for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
356       {
357         Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
358         //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
359         //Distance to bad module 1
360         snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
361         snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
362                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
363         fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
364         fhRe1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
365         fhRe1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
366         //printf("name: %s\n ",fhRe1[index]->GetName());
367         outputContainer->Add(fhRe1[index]) ;
368         
369         if(fFillBadDistHisto)
370         {
371           //Distance to bad module 2
372           snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
373           snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
374                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
375           fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
376           fhRe2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
377           fhRe2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
378           outputContainer->Add(fhRe2[index]) ;
379           
380           //Distance to bad module 3
381           snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
382           snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
383                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
384           fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
385           fhRe3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
386           fhRe3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
387           outputContainer->Add(fhRe3[index]) ;
388         }
389         
390         //Inverse pT
391         if(fMakeInvPtPlots)
392         {
393           //Distance to bad module 1
394           snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
395           snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
396                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
397           fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
398           fhReInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
399           fhReInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
400           outputContainer->Add(fhReInvPt1[index]) ;
401           
402           if(fFillBadDistHisto){
403             //Distance to bad module 2
404             snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
405             snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
406                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
407             fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
408             fhReInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
409             fhReInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
410             outputContainer->Add(fhReInvPt2[index]) ;
411             
412             //Distance to bad module 3
413             snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
414             snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
415                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
416             fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
417             fhReInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
418             fhReInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
419             outputContainer->Add(fhReInvPt3[index]) ;
420           }
421         }
422         
423         if(DoOwnMix())
424         {
425           //Distance to bad module 1
426           snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
427           snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
428                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
429           fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
430           fhMi1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
431           fhMi1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
432           outputContainer->Add(fhMi1[index]) ;
433           if(fFillBadDistHisto){
434             //Distance to bad module 2
435             snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
436             snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
437                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
438             fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
439             fhMi2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
440             fhMi2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
441             outputContainer->Add(fhMi2[index]) ;
442             
443             //Distance to bad module 3
444             snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
445             snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
446                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
447             fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
448             fhMi3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
449             fhMi3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
450             outputContainer->Add(fhMi3[index]) ;
451           }
452           
453           //Inverse pT
454           if(fMakeInvPtPlots)
455           {
456             //Distance to bad module 1
457             snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
458             snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
459                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
460             fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
461             fhMiInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
462             fhMiInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
463             outputContainer->Add(fhMiInvPt1[index]) ;
464             if(fFillBadDistHisto){
465               //Distance to bad module 2
466               snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
467               snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
468                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
469               fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
470               fhMiInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
471               fhMiInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
472               outputContainer->Add(fhMiInvPt2[index]) ;
473               
474               //Distance to bad module 3
475               snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
476               snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
477                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
478               fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
479               fhMiInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
480               fhMiInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
481               outputContainer->Add(fhMiInvPt3[index]) ;
482             }
483           }
484         }
485       }
486     }
487   }
488   
489   if(fFillAsymmetryHisto)
490   {
491     fhRePtAsym = new TH2F("hRePtAsym","#it{Asymmetry} vs #it{p}_{T} , for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
492     fhRePtAsym->SetXTitle("#it{p}_{T} (GeV/#it{c})");
493     fhRePtAsym->SetYTitle("#it{Asymmetry}");
494     outputContainer->Add(fhRePtAsym);
495     
496     fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
497     fhRePtAsymPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
498     fhRePtAsymPi0->SetYTitle("Asymmetry");
499     outputContainer->Add(fhRePtAsymPi0);
500     
501     fhRePtAsymEta = new TH2F("hRePtAsymEta","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
502     fhRePtAsymEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
503     fhRePtAsymEta->SetYTitle("#it{Asymmetry}");
504     outputContainer->Add(fhRePtAsymEta);
505   }
506   
507   if(fMultiCutAna)
508   {
509     fhRePIDBits         = new TH2F*[fNPIDBits];
510     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
511       snprintf(key,   buffersize,"hRe_pidbit%d",ipid) ;
512       snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
513       fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
514       fhRePIDBits[ipid]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
515       fhRePIDBits[ipid]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
516       outputContainer->Add(fhRePIDBits[ipid]) ;
517     }// pid bit loop
518     
519     fhRePtNCellAsymCuts    = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
520     fhMiPtNCellAsymCuts    = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
521     
522     if(fFillSMCombinations)
523       for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
524     
525     
526     for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
527     {
528       for(Int_t icell=0; icell<fNCellNCuts; icell++)
529       {
530         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
531         {
532           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
533           snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
534           Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
535           //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
536           fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
537           fhRePtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
538           fhRePtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
539           outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
540           
541           snprintf(key,   buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
542           snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
543           fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
544           fhMiPtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
545           fhMiPtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
546           outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
547           
548           if(fFillSMCombinations)
549           {
550             for(Int_t iSM = 0; iSM < fNModules; iSM++)
551             {
552               snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
553               snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
554                        fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
555               fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
556               fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
557               fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
558               outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
559             }
560             
561           }
562         }
563       }
564     }
565     
566     if(ntrmbins!=0)
567     {
568       fhRePtMult = new TH3F*[fNAsymCuts] ;
569       for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
570       {
571         fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(#it{p}_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
572                                      nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
573         fhRePtMult[iasym]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
574         fhRePtMult[iasym]->SetYTitle("Track multiplicity");
575         fhRePtMult[iasym]->SetZTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
576         outputContainer->Add(fhRePtMult[iasym]) ;
577       }
578     }
579   }// multi cuts analysis
580   
581   if(fFillSSCombinations)
582   {
583     
584     fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
585                          nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
586     fhReSS[0]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
587     fhReSS[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
588     outputContainer->Add(fhReSS[0]) ;
589     
590     
591     fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
592                          nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
593     fhReSS[1]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
594     fhReSS[1]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
595     outputContainer->Add(fhReSS[1]) ;
596     
597     
598     fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
599                          nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
600     fhReSS[2]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
601     fhReSS[2]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
602     outputContainer->Add(fhReSS[2]) ;
603   }
604   
605   if(DoOwnMix())
606   {
607     fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
608                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
609                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
610     fhEventBin->SetXTitle("bin");
611     outputContainer->Add(fhEventBin) ;
612     
613     
614     fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
615                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
616                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
617     fhEventMixBin->SetXTitle("bin");
618     outputContainer->Add(fhEventMixBin) ;
619   }
620   
621   if(GetNCentrBin()>1)
622   {
623     fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
624     fhCentrality->SetXTitle("Centrality bin");
625     outputContainer->Add(fhCentrality) ;
626     
627     fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
628     fhCentralityNoPair->SetXTitle("Centrality bin");
629     outputContainer->Add(fhCentralityNoPair) ;
630   }
631   
632   if(GetNRPBin() > 1 && GetNCentrBin()>1 )
633   {
634     fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
635     fhEventPlaneResolution->SetYTitle("Resolution");
636     fhEventPlaneResolution->SetXTitle("Centrality Bin");
637     outputContainer->Add(fhEventPlaneResolution) ;
638   }
639   
640   if(fFillAngleHisto)
641   {
642     fhRealOpeningAngle  = new TH2F
643     ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
644     fhRealOpeningAngle->SetYTitle("#theta(rad)");
645     fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
646     outputContainer->Add(fhRealOpeningAngle) ;
647     
648     fhRealCosOpeningAngle  = new TH2F
649     ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
650     fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
651     fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
652     outputContainer->Add(fhRealCosOpeningAngle) ;
653     
654     if(DoOwnMix())
655     {
656       fhMixedOpeningAngle  = new TH2F
657       ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
658       fhMixedOpeningAngle->SetYTitle("#theta(rad)");
659       fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
660       outputContainer->Add(fhMixedOpeningAngle) ;
661       
662       fhMixedCosOpeningAngle  = new TH2F
663       ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
664       fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
665       fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
666       outputContainer->Add(fhMixedCosOpeningAngle) ;
667       
668     }
669   }
670   
671   //Histograms filled only if MC data is requested
672   if(IsDataMC())
673   {
674     fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
675                                     nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
676     fhReMCFromConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
677     fhReMCFromConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
678     outputContainer->Add(fhReMCFromConversion) ;
679     
680     fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
681                                        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
682     fhReMCFromNotConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
683     fhReMCFromNotConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
684     outputContainer->Add(fhReMCFromNotConversion) ;
685     
686     fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
687                                        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
688     fhReMCFromMixConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
689     fhReMCFromMixConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
690     outputContainer->Add(fhReMCFromMixConversion) ;
691     
692     //Pi0
693     
694     fhPrimPi0E     = new TH1F("hPrimPi0E","Primary #pi^{0} E, |#it{Y}|<1",nptbins,ptmin,ptmax) ;
695     fhPrimPi0AccE  = new TH1F("hPrimPi0AccE","Primary #pi^{0} #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
696     fhPrimPi0E   ->SetXTitle("#it{E} (GeV)");
697     fhPrimPi0AccE->SetXTitle("#it{E} (GeV)");
698     outputContainer->Add(fhPrimPi0E) ;
699     outputContainer->Add(fhPrimPi0AccE) ;
700     
701     fhPrimPi0Pt     = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , |#it{Y}|<1",nptbins,ptmin,ptmax) ;
702     fhPrimPi0AccPt  = new TH1F("hPrimPi0AccPt","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
703     fhPrimPi0Pt   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
704     fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
705     outputContainer->Add(fhPrimPi0Pt) ;
706     outputContainer->Add(fhPrimPi0AccPt) ;
707     
708     Int_t netabinsopen =  TMath::Nint(netabins*4/(etamax-etamin));
709     fhPrimPi0Y      = new TH2F("hPrimPi0Rapidity","Rapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
710     fhPrimPi0Y   ->SetYTitle("#it{Rapidity}");
711     fhPrimPi0Y   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
712     outputContainer->Add(fhPrimPi0Y) ;
713     
714     fhPrimPi0Yeta      = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
715     fhPrimPi0Yeta   ->SetYTitle("#eta");
716     fhPrimPi0Yeta   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
717     outputContainer->Add(fhPrimPi0Yeta) ;
718     
719     fhPrimPi0YetaYcut      = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
720     fhPrimPi0YetaYcut   ->SetYTitle("#eta");
721     fhPrimPi0YetaYcut   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
722     outputContainer->Add(fhPrimPi0YetaYcut) ;
723     
724     fhPrimPi0AccY   = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
725     fhPrimPi0AccY->SetYTitle("Rapidity");
726     fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
727     outputContainer->Add(fhPrimPi0AccY) ;
728     
729     fhPrimPi0AccYeta      = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
730     fhPrimPi0AccYeta   ->SetYTitle("#eta");
731     fhPrimPi0AccYeta   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
732     outputContainer->Add(fhPrimPi0AccYeta) ;
733     
734     Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
735     fhPrimPi0Phi    = new TH2F("hPrimPi0Phi","#phi of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
736     fhPrimPi0Phi->SetYTitle("#phi (deg)");
737     fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
738     outputContainer->Add(fhPrimPi0Phi) ;
739     
740     fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","#phi of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,
741                                nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
742     fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
743     fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
744     outputContainer->Add(fhPrimPi0AccPhi) ;
745     
746     fhPrimPi0PtCentrality     = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, |#it{Y}|<1",
747                                          nptbins,ptmin,ptmax, 100, 0, 100) ;
748     fhPrimPi0AccPtCentrality  = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality",
749                                          nptbins,ptmin,ptmax, 100, 0, 100) ;
750     fhPrimPi0PtCentrality   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
751     fhPrimPi0AccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
752     fhPrimPi0PtCentrality   ->SetYTitle("Centrality");
753     fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
754     outputContainer->Add(fhPrimPi0PtCentrality) ;
755     outputContainer->Add(fhPrimPi0AccPtCentrality) ;
756     
757     fhPrimPi0PtEventPlane     = new TH2F("hPrimPi0PtEventPlane","Primary #pi^{0} #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",
758                                          nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
759     fhPrimPi0AccPtEventPlane  = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle",
760                                          nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
761     fhPrimPi0PtEventPlane   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
762     fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
763     fhPrimPi0PtEventPlane   ->SetYTitle("Event Plane Angle (rad)");
764     fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
765     outputContainer->Add(fhPrimPi0PtEventPlane) ;
766     outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
767     
768     //Eta
769     
770     fhPrimEtaE     = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
771     fhPrimEtaAccE  = new TH1F("hPrimEtaAccE","Primary #eta #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
772     fhPrimEtaE   ->SetXTitle("#it{E} (GeV)");
773     fhPrimEtaAccE->SetXTitle("#it{E} (GeV)");
774     outputContainer->Add(fhPrimEtaE) ;
775     outputContainer->Add(fhPrimEtaAccE) ;
776     
777     fhPrimEtaPt     = new TH1F("hPrimEtaPt","Primary #eta #it{p}_{T}",nptbins,ptmin,ptmax) ;
778     fhPrimEtaAccPt  = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
779     fhPrimEtaPt   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
780     fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
781     outputContainer->Add(fhPrimEtaPt) ;
782     outputContainer->Add(fhPrimEtaAccPt) ;
783     
784     fhPrimEtaY      = new TH2F("hPrimEtaRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
785     fhPrimEtaY->SetYTitle("#it{Rapidity}");
786     fhPrimEtaY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
787     outputContainer->Add(fhPrimEtaY) ;
788     
789     fhPrimEtaYeta      = new TH2F("hPrimEtaPseudoRapidityEta","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
790     fhPrimEtaYeta->SetYTitle("#it{Rapidity}");
791     fhPrimEtaYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
792     outputContainer->Add(fhPrimEtaYeta) ;
793     
794     fhPrimEtaYetaYcut      = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary #eta, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
795     fhPrimEtaYetaYcut->SetYTitle("#it{Pseudorapidity}");
796     fhPrimEtaYetaYcut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
797     outputContainer->Add(fhPrimEtaYetaYcut) ;
798     
799     fhPrimEtaAccY   = new TH2F("hPrimEtaAccRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
800     fhPrimEtaAccY->SetYTitle("#it{Rapidity}");
801     fhPrimEtaAccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
802     outputContainer->Add(fhPrimEtaAccY) ;
803     
804     fhPrimEtaAccYeta  = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
805     fhPrimEtaAccYeta->SetYTitle("#it{Pseudorapidity}");
806     fhPrimEtaAccYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
807     outputContainer->Add(fhPrimEtaAccYeta) ;
808     
809     fhPrimEtaPhi    = new TH2F("hPrimEtaPhi","Azimuthal of primary #eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
810     fhPrimEtaPhi->SetYTitle("#phi (deg)");
811     fhPrimEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
812     outputContainer->Add(fhPrimEtaPhi) ;
813     
814     fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary #eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
815     fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
816     fhPrimEtaAccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
817     outputContainer->Add(fhPrimEtaAccPhi) ;
818     
819     fhPrimEtaPtCentrality     = new TH2F("hPrimEtaPtCentrality","Primary #eta #it{p}_{T} vs reco centrality, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
820     fhPrimEtaAccPtCentrality  = new TH2F("hPrimEtaAccPtCentrality","Primary #eta with both photons in acceptance #it{p}_{T} vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
821     fhPrimEtaPtCentrality   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
822     fhPrimEtaAccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
823     fhPrimEtaPtCentrality   ->SetYTitle("Centrality");
824     fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
825     outputContainer->Add(fhPrimEtaPtCentrality) ;
826     outputContainer->Add(fhPrimEtaAccPtCentrality) ;
827     
828     fhPrimEtaPtEventPlane     = new TH2F("hPrimEtaPtEventPlane","Primary #eta #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
829     fhPrimEtaAccPtEventPlane  = new TH2F("hPrimEtaAccPtEventPlane","Primary #eta with both #gamma_{decay} in acceptance #it{p}_{T} vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
830     fhPrimEtaPtEventPlane   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
831     fhPrimEtaAccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
832     fhPrimEtaPtEventPlane   ->SetYTitle("Event Plane Angle (rad)");
833     fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
834     outputContainer->Add(fhPrimEtaPtEventPlane) ;
835     outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
836     
837     if(fFillAngleHisto)
838     {
839       fhPrimPi0OpeningAngle  = new TH2F
840       ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
841       fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
842       fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
843       outputContainer->Add(fhPrimPi0OpeningAngle) ;
844       
845       fhPrimPi0OpeningAngleAsym  = new TH2F
846       ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5);
847       fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
848       fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
849       outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
850       
851       fhPrimPi0CosOpeningAngle  = new TH2F
852       ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
853       fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
854       fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
855       outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
856       
857       fhPrimEtaOpeningAngle  = new TH2F
858       ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,0,0.5);
859       fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
860       fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
861       outputContainer->Add(fhPrimEtaOpeningAngle) ;
862       
863       fhPrimEtaOpeningAngleAsym  = new TH2F
864       ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5);
865       fhPrimEtaOpeningAngleAsym->SetXTitle("|#it{A}|=| (#it{E}_{1}-#it{E}_{2}) / (#it{E}_{1}+#it{E}_{2}) |");
866       fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
867       outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
868       
869       
870       fhPrimEtaCosOpeningAngle  = new TH2F
871       ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
872       fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
873       fhPrimEtaCosOpeningAngle->SetXTitle("#it{E}_{ #eta} (GeV)");
874       outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
875     }
876     
877     if(fFillOriginHisto)
878     {
879       //Prim origin
880       //Pi0
881       fhPrimPi0PtOrigin     = new TH2F("hPrimPi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
882       fhPrimPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
883       fhPrimPi0PtOrigin->SetYTitle("Origin");
884       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
885       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
886       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
887       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
888       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
889       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
890       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
891       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
892       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
893       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
894       outputContainer->Add(fhPrimPi0PtOrigin) ;
895       
896       fhMCPi0PtOrigin     = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
897       fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
898       fhMCPi0PtOrigin->SetYTitle("Origin");
899       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
900       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
901       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
902       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
903       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
904       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
905       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
906       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
907       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
908       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
909       outputContainer->Add(fhMCPi0PtOrigin) ;
910       
911       //Eta
912       fhPrimEtaPtOrigin     = new TH2F("hPrimEtaPtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
913       fhPrimEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
914       fhPrimEtaPtOrigin->SetYTitle("Origin");
915       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
916       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
917       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
918       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
919       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
920       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
921       
922       outputContainer->Add(fhPrimEtaPtOrigin) ;
923       
924       fhMCEtaPtOrigin     = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
925       fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
926       fhMCEtaPtOrigin->SetYTitle("Origin");
927       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
928       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
929       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
930       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
931       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
932       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
933       
934       outputContainer->Add(fhMCEtaPtOrigin) ;
935       
936       fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",
937                                    200,0.,20.,5000,0,500) ;
938       fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
939       fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
940       outputContainer->Add(fhMCPi0ProdVertex) ;
941       
942       fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",
943                                    200,0.,20.,5000,0,500) ;
944       fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
945       fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
946       outputContainer->Add(fhMCEtaProdVertex) ;
947       
948       fhPrimPi0ProdVertex = new TH2F("hPrimPi0ProdVertex","generated #pi^{0} #it{p}_{T} vs production vertex",
949                                      200,0.,20.,5000,0,500) ;
950       fhPrimPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
951       fhPrimPi0ProdVertex->SetYTitle("#it{R} (cm)");
952       outputContainer->Add(fhPrimPi0ProdVertex) ;
953       
954       fhPrimEtaProdVertex = new TH2F("hPrimEtaProdVertex","generated #eta #it{p}_{T} vs production vertex",
955                                      200,0.,20.,5000,0,500) ;
956       fhPrimEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
957       fhPrimEtaProdVertex->SetYTitle("#it{R} (cm)");
958       outputContainer->Add(fhPrimEtaProdVertex) ;
959       
960       for(Int_t i = 0; i<13; i++)
961       {
962         fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
963         fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
964         fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
965         outputContainer->Add(fhMCOrgMass[i]) ;
966         
967         fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
968         fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
969         fhMCOrgAsym[i]->SetYTitle("A");
970         outputContainer->Add(fhMCOrgAsym[i]) ;
971         
972         fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
973         fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
974         fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
975         outputContainer->Add(fhMCOrgDeltaEta[i]) ;
976         
977         fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
978         fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
979         fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
980         outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
981         
982       }
983       
984       if(fMultiCutAnaSim)
985       {
986         fhMCPi0MassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
987         fhMCPi0MassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
988         fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
989         fhMCEtaMassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
990         fhMCEtaMassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
991         fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
992         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
993           for(Int_t icell=0; icell<fNCellNCuts; icell++){
994             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
995               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
996               
997               fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
998                                                  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
999                                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1000               fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1001               fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1002               outputContainer->Add(fhMCPi0MassPtRec[index]) ;
1003               
1004               fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1005                                                   Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1006                                                   nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1007               fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1008               fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1009               outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
1010               
1011               fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1012                                                    Form("Generated vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1013                                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1014               fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1015               fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1016               outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
1017               
1018               fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1019                                                  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1020                                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1021               fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1022               fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1023               outputContainer->Add(fhMCEtaMassPtRec[index]) ;
1024               
1025               fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1026                                                   Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1027                                                   nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1028               fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1029               fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1030               outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1031               
1032               fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1033                                                    Form("Generated vs reconstructed #it{p}_{T} of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1034                                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1035               fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1036               fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1037               outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1038             }
1039           }
1040         }
1041       }//multi cut ana
1042       else
1043       {
1044         fhMCPi0MassPtTrue  = new TH2F*[1];
1045         fhMCPi0PtTruePtRec = new TH2F*[1];
1046         fhMCEtaMassPtTrue  = new TH2F*[1];
1047         fhMCEtaPtTruePtRec = new TH2F*[1];
1048         
1049         fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1050         fhMCPi0MassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1051         fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1052         outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1053         
1054         fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1055         fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1056         fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1057         outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1058         
1059         fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1060         fhMCEtaMassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1061         fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1062         outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1063         
1064         fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1065         fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1066         fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1067         outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1068       }
1069     }
1070   }
1071   
1072   if(fFillSMCombinations)
1073   {
1074     TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1075     for(Int_t imod=0; imod<fNModules; imod++)
1076     {
1077       //Module dependent invariant mass
1078       snprintf(key, buffersize,"hReMod_%d",imod) ;
1079       snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1080       fhReMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1081       fhReMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1082       fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1083       outputContainer->Add(fhReMod[imod]) ;
1084       if(GetCalorimeter()==kPHOS)
1085       {
1086         snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1087         snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1088         fhReDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1089         fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1090         fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1091         outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1092       }
1093       else
1094       {//EMCAL
1095         if(imod<fNModules/2)
1096         {
1097           snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1098           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1099           fhReSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1100           fhReSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1101           fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1102           outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1103         }
1104         if(imod<fNModules-2)
1105         {
1106           snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1107           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1108           fhReSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1109           fhReSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1110           fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1111           outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1112         }
1113       }//EMCAL
1114       
1115       if(DoOwnMix())
1116       {
1117         snprintf(key, buffersize,"hMiMod_%d",imod) ;
1118         snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1119         fhMiMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1120         fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1121         fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1122         outputContainer->Add(fhMiMod[imod]) ;
1123         
1124         if(GetCalorimeter()==kPHOS)
1125         {
1126           snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1127           snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1128           fhMiDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1129           fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1130           fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1131           outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1132         }//PHOS
1133         else
1134         {//EMCAL
1135           if(imod<fNModules/2)
1136           {
1137             snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1138             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1139             fhMiSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1140             fhMiSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1141             fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1142             outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1143           }
1144           if(imod<fNModules-2){
1145             
1146             snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1147             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1148             fhMiSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1149             fhMiSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1150             fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1151             outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1152           }
1153         }//EMCAL
1154       }// own mix
1155     }//loop combinations
1156   } // SM combinations
1157   
1158   if(fFillArmenterosThetaStar && IsDataMC())
1159   {
1160     TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
1161     Int_t narmbins = 400;
1162     Float_t armmin = 0;
1163     Float_t armmax = 0.4;
1164     
1165     for(Int_t i = 0; i < 4; i++)
1166     {
1167       fhArmPrimPi0[i] =  new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
1168                                   Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
1169                                   200, -1, 1, narmbins,armmin,armmax);
1170       fhArmPrimPi0[i]->SetYTitle("#it{p}_{T}^{Arm}");
1171       fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
1172       outputContainer->Add(fhArmPrimPi0[i]) ;
1173       
1174       fhArmPrimEta[i] =  new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
1175                                   Form("Armenteros of primary #eta, %s",ebin[i].Data()),
1176                                   200, -1, 1, narmbins,armmin,armmax);
1177       fhArmPrimEta[i]->SetYTitle("#it{p}_{T}^{Arm}");
1178       fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
1179       outputContainer->Add(fhArmPrimEta[i]) ;
1180       
1181     }
1182     
1183     // Same as asymmetry ...
1184     fhCosThStarPrimPi0  = new TH2F
1185     ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
1186     fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
1187     fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
1188     outputContainer->Add(fhCosThStarPrimPi0) ;
1189     
1190     fhCosThStarPrimEta  = new TH2F
1191     ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
1192     fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
1193     fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
1194     outputContainer->Add(fhCosThStarPrimEta) ;
1195     
1196   }
1197   
1198   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1199   //
1200   //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1201   //
1202   //  }
1203   
1204   return outputContainer;
1205 }
1206
1207 //___________________________________________________
1208 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1209 {
1210   //Print some relevant parameters set for the analysis
1211   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1212   AliAnaCaloTrackCorrBaseClass::Print(" ");
1213   
1214   printf("Number of bins in Centrality:  %d \n",GetNCentrBin()) ;
1215   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1216   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1217   printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1218   printf("Pair in same Module: %d \n",fSameSM) ;
1219   printf("Cuts: \n") ;
1220   // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1221   printf("Number of modules:             %d \n",fNModules) ;
1222   printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1223   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1224   printf("\tasymmetry < ");
1225   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1226   printf("\n");
1227   
1228   printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1229   printf("\tPID bit = ");
1230   for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1231   printf("\n");
1232   
1233   if(fMultiCutAna)
1234   {
1235     printf("pT cuts: n = %d, \n",fNPtCuts) ;
1236     printf("\tpT > ");
1237     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1238     printf("GeV/c\n");
1239     
1240     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1241     printf("\tnCell > ");
1242     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1243     printf("\n");
1244     
1245   }
1246   printf("------------------------------------------------------\n") ;
1247 }
1248
1249 //________________________________________
1250 void AliAnaPi0::FillAcceptanceHistograms()
1251 {
1252   //Fill acceptance histograms if MC data is available
1253   
1254   Double_t mesonY   = -100 ;
1255   Double_t mesonE   = -1 ;
1256   Double_t mesonPt  = -1 ;
1257   Double_t mesonPhi =  100 ;
1258   Double_t mesonYeta= -1 ;
1259   
1260   Int_t    pdg     = 0 ;
1261   Int_t    nprim   = 0 ;
1262   Int_t    nDaught = 0 ;
1263   Int_t    iphot1  = 0 ;
1264   Int_t    iphot2  = 0 ;
1265   
1266   Float_t cen = GetEventCentrality();
1267   Float_t ep  = GetEventPlaneAngle();
1268   
1269   TParticle        * primStack = 0;
1270   AliAODMCParticle * primAOD   = 0;
1271   
1272   // Get the ESD MC particles container
1273   AliStack * stack = 0;
1274   if( GetReader()->ReadStack() )
1275   {
1276     stack = GetMCStack();
1277     if(!stack ) return;
1278     nprim = stack->GetNtrack();
1279   }
1280   
1281   // Get the AOD MC particles container
1282   TClonesArray * mcparticles = 0;
1283   if( GetReader()->ReadAODMCParticles() )
1284   {
1285     mcparticles = GetReader()->GetAODMCParticles();
1286     if( !mcparticles ) return;
1287     nprim = mcparticles->GetEntriesFast();
1288   }
1289   
1290   for(Int_t i=0 ; i < nprim; i++)
1291   {
1292     if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1293     
1294     if(GetReader()->ReadStack())
1295     {
1296       primStack = stack->Particle(i) ;
1297       if(!primStack)
1298       {
1299         AliWarning("ESD primaries pointer not available!!");
1300         continue;
1301       }
1302       
1303       // If too small  skip
1304       if( primStack->Energy() < 0.4 ) continue;
1305       
1306       pdg       = primStack->GetPdgCode();
1307       nDaught   = primStack->GetNDaughters();
1308       iphot1    = primStack->GetDaughter(0) ;
1309       iphot2    = primStack->GetDaughter(1) ;
1310       if(primStack->Energy() == TMath::Abs(primStack->Pz()))  continue ; //Protection against floating point exception
1311       
1312       //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1313       //       prim->GetName(), prim->GetPdgCode());
1314       
1315       //Photon kinematics
1316       primStack->Momentum(fPi0Mom);
1317       
1318       mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
1319     }
1320     else
1321     {
1322       primAOD = (AliAODMCParticle *) mcparticles->At(i);
1323       if(!primAOD)
1324       {
1325         AliWarning("AOD primaries pointer not available!!");
1326         continue;
1327       }
1328       
1329       // If too small  skip
1330       if( primAOD->E() < 0.4 ) continue;
1331       
1332       pdg     = primAOD->GetPdgCode();
1333       nDaught = primAOD->GetNDaughters();
1334       iphot1  = primAOD->GetFirstDaughter() ;
1335       iphot2  = primAOD->GetLastDaughter() ;
1336       
1337       if(primAOD->E() == TMath::Abs(primAOD->Pz()))  continue ; //Protection against floating point exception
1338       
1339       //Photon kinematics
1340       fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1341       
1342       mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
1343     }
1344     
1345     // Select only pi0 or eta
1346     if( pdg != 111 && pdg != 221) continue ;
1347     
1348     mesonPt  = fPi0Mom.Pt () ;
1349     mesonE   = fPi0Mom.E  () ;
1350     mesonYeta= fPi0Mom.Eta() ;
1351     mesonPhi = fPi0Mom.Phi() ;
1352     if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
1353     mesonPhi *= TMath::RadToDeg();
1354     
1355     if(pdg == 111)
1356     {
1357       if(TMath::Abs(mesonY) < 1.0)
1358       {
1359         fhPrimPi0E  ->Fill(mesonE ) ;
1360         fhPrimPi0Pt ->Fill(mesonPt) ;
1361         fhPrimPi0Phi->Fill(mesonPt, mesonPhi) ;
1362         
1363         fhPrimPi0YetaYcut    ->Fill(mesonPt,mesonYeta) ;
1364         fhPrimPi0PtCentrality->Fill(mesonPt,cen) ;
1365         fhPrimPi0PtEventPlane->Fill(mesonPt,ep ) ;
1366       }
1367       
1368       fhPrimPi0Y   ->Fill(mesonPt, mesonY) ;
1369       fhPrimPi0Yeta->Fill(mesonPt, mesonYeta) ;
1370     }
1371     else if(pdg == 221)
1372     {
1373       if(TMath::Abs(mesonY) < 1.0)
1374       {
1375         fhPrimEtaE  ->Fill(mesonE ) ;
1376         fhPrimEtaPt ->Fill(mesonPt) ;
1377         fhPrimEtaPhi->Fill(mesonPt, mesonPhi) ;
1378         
1379         fhPrimEtaYetaYcut    ->Fill(mesonPt,mesonYeta) ;
1380         fhPrimEtaPtCentrality->Fill(mesonPt,cen) ;
1381         fhPrimEtaPtEventPlane->Fill(mesonPt,ep ) ;
1382       }
1383       
1384       fhPrimEtaY   ->Fill(mesonPt, mesonY) ;
1385       fhPrimEtaYeta->Fill(mesonPt, mesonYeta) ;
1386     }
1387     
1388     //Origin of meson
1389     if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
1390     {
1391       Int_t momindex  = -1;
1392       Int_t mompdg    = -1;
1393       Int_t momstatus = -1;
1394       Float_t momR    =  0;
1395       if(GetReader()->ReadStack())          momindex = primStack->GetFirstMother();
1396       if(GetReader()->ReadAODMCParticles()) momindex = primAOD  ->GetMother();
1397       
1398       if(momindex >= 0 && momindex < nprim)
1399       {
1400         if(GetReader()->ReadStack())
1401         {
1402           TParticle* mother = stack->Particle(momindex);
1403           mompdg    = TMath::Abs(mother->GetPdgCode());
1404           momstatus = mother->GetStatusCode();
1405           momR      = mother->R();
1406         }
1407         
1408         if(GetReader()->ReadAODMCParticles())
1409         {
1410           AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex);
1411           mompdg    = TMath::Abs(mother->GetPdgCode());
1412           momstatus = mother->GetStatus();
1413           momR      = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1414         }
1415         
1416         if(pdg == 111)
1417         {
1418           if     (momstatus  == 21)fhPrimPi0PtOrigin->Fill(mesonPt,0.5);//parton
1419           else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt,1.5);//quark
1420           else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(mesonPt,2.5);// resonances
1421           else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(mesonPt,8.5);//eta
1422           else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(mesonPt,9.5);//eta prime
1423           else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(mesonPt,4.5);//rho
1424           else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(mesonPt,5.5);//omega
1425           else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0S, k+-,k*
1426           else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0L
1427           else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(mesonPt,3.5);//resonances
1428           else                      fhPrimPi0PtOrigin->Fill(mesonPt,7.5);//other?
1429           
1430           //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1431           //                   momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
1432           
1433           fhPrimPi0ProdVertex->Fill(mesonPt,momR);
1434           
1435         }//pi0
1436         else
1437         {
1438           if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt,0.5);//parton
1439           else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(mesonPt,1.5);//quark
1440           else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(mesonPt,2.5);//qq resonances
1441           else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(mesonPt,5.5);//eta prime
1442           else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(mesonPt,3.5);//resonances
1443           else fhPrimEtaPtOrigin->Fill(mesonPt,4.5);//stable, conversions?
1444           //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1445           
1446           fhPrimEtaProdVertex->Fill(mesonPt,momR);
1447           
1448         }
1449       } // pi0 has mother
1450     }
1451     
1452     //Check if both photons hit Calorimeter
1453     if(nDaught != 2 ) continue; //Only interested in 2 gamma decay
1454     
1455     if(iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim) continue ;
1456     
1457     Int_t pdg1 = 0;
1458     Int_t pdg2 = 0;
1459     Bool_t inacceptance1 = kTRUE;
1460     Bool_t inacceptance2 = kTRUE;
1461     
1462     if(GetReader()->ReadStack())
1463     {
1464       TParticle * phot1 = stack->Particle(iphot1) ;
1465       TParticle * phot2 = stack->Particle(iphot2) ;
1466       
1467       if(!phot1 || !phot2) continue ;
1468       
1469       pdg1 = phot1->GetPdgCode();
1470       pdg2 = phot2->GetPdgCode();
1471       
1472       phot1->Momentum(fPhotonMom1);
1473       phot2->Momentum(fPhotonMom2);
1474       
1475       // Check if photons hit the Calorimeter acceptance
1476       if(IsRealCaloAcceptanceOn())
1477       {
1478         if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1479         if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1480       }
1481     }
1482     
1483     if(GetReader()->ReadAODMCParticles())
1484     {
1485       AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ;
1486       AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ;
1487       
1488       if(!phot1 || !phot2) continue ;
1489       
1490       pdg1 = phot1->GetPdgCode();
1491       pdg2 = phot2->GetPdgCode();
1492       
1493       fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1494       fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1495       
1496       // Check if photons hit the Calorimeter acceptance
1497       if(IsRealCaloAcceptanceOn())
1498       {
1499         if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1500         if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1501       }
1502     }
1503     
1504     if( pdg1 != 22 || pdg2 !=22) continue ;
1505     
1506     // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
1507     if(IsFiducialCutOn())
1508     {
1509       if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
1510       if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
1511     }
1512     
1513     if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg);
1514     
1515     if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
1516     {
1517       Int_t absID1=0;
1518       Int_t absID2=0;
1519       
1520       Float_t photonPhi1 = fPhotonMom1.Phi();
1521       Float_t photonPhi2 = fPhotonMom2.Phi();
1522       
1523       if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
1524       if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
1525       
1526       GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
1527       GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
1528       
1529       Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1530       Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1531       
1532       Int_t j=0;
1533       Bool_t sameSector = kFALSE;
1534       for(Int_t isector = 0; isector < fNModules/2; isector++)
1535       {
1536         j=2*isector;
1537         if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1538       }
1539       
1540       if(sm1!=sm2 && !sameSector)
1541       {
1542         inacceptance1 = kFALSE;
1543         inacceptance2 = kFALSE;
1544       }
1545       //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1546       //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1547       //                    inacceptance = kTRUE;
1548     }
1549     
1550     AliDebug(2,Form("Accepted in %s?: m (%2.2f,%2.2f,%2.2f), p1 (%2.2f,%2.2f,%2.2f), p2 (%2.2f,%2.2f,%2.2f) : in1 %d, in2 %d",
1551                     GetCalorimeterString().Data(),
1552                     mesonPt,mesonYeta,mesonPhi,
1553                     fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
1554                     fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
1555                     inacceptance1, inacceptance2));
1556     
1557     if(inacceptance1 && inacceptance2)
1558     {
1559       Float_t  asym  = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
1560       Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
1561       
1562       AliDebug(2,Form("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f",pdg,mesonPt,mesonPhi,mesonYeta));
1563       
1564       if(pdg==111)
1565       {
1566         fhPrimPi0AccE   ->Fill(mesonE) ;
1567         fhPrimPi0AccPt  ->Fill(mesonPt) ;
1568         fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi) ;
1569         fhPrimPi0AccY   ->Fill(mesonPt, mesonY) ;
1570         fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta) ;
1571         fhPrimPi0AccPtCentrality->Fill(mesonPt,cen) ;
1572         fhPrimPi0AccPtEventPlane->Fill(mesonPt,ep ) ;
1573         
1574         if(fFillAngleHisto)
1575         {
1576           fhPrimPi0OpeningAngle    ->Fill(mesonPt,angle);
1577           if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1578           fhPrimPi0CosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle));
1579         }
1580       }
1581       else if(pdg==221)
1582       {
1583         fhPrimEtaAccE   ->Fill(mesonE ) ;
1584         fhPrimEtaAccPt  ->Fill(mesonPt) ;
1585         fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi) ;
1586         fhPrimEtaAccY   ->Fill(mesonPt, mesonY) ;
1587         fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta) ;
1588         fhPrimEtaAccPtCentrality->Fill(mesonPt,cen) ;
1589         fhPrimEtaAccPtEventPlane->Fill(mesonPt,ep ) ;
1590         
1591         if(fFillAngleHisto)
1592         {
1593           fhPrimEtaOpeningAngle    ->Fill(mesonPt,angle);
1594           if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1595           fhPrimEtaCosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle));
1596         }
1597       }
1598     }//Accepted
1599     
1600   }//loop on primaries
1601   
1602 }
1603
1604 //________________________________________________
1605 void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg)
1606 {
1607   // Fill armenteros plots
1608   
1609   // Get pTArm and AlphaArm
1610   Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P();
1611   Float_t momentumDaughter1AlongMother = 0.;
1612   Float_t momentumDaughter2AlongMother = 0.;
1613   
1614   if (momentumSquaredMother > 0.)
1615   {
1616     momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1617     momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1618   }
1619   
1620   Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
1621   Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1622   
1623   Float_t pTArm = 0.;
1624   if (ptArmSquared > 0.)
1625     pTArm = sqrt(ptArmSquared);
1626   
1627   Float_t alphaArm = 0.;
1628   if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
1629     alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1630   
1631   fPhotonMom1Boost = fPhotonMom1;
1632   fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector());
1633   Float_t  cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect()));
1634   
1635   Float_t en   = fPi0Mom.Energy();
1636   Int_t   ebin = -1;
1637   if(en > 8  && en <= 12) ebin = 0;
1638   if(en > 12 && en <= 16) ebin = 1;
1639   if(en > 16 && en <= 20) ebin = 2;
1640   if(en > 20)             ebin = 3;
1641   if(ebin < 0 || ebin > 3) return ;
1642   
1643   
1644   if(pdg==111)
1645   {
1646     fhCosThStarPrimPi0->Fill(en,cosThStar);
1647     fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
1648   }
1649   else
1650   {
1651     fhCosThStarPrimEta->Fill(en,cosThStar);
1652     fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
1653   }
1654   
1655   //  if(GetDebug() > 2 )
1656   //  {
1657   //    Float_t asym = 0;
1658   //    if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
1659   //
1660   //    printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1661   //         en,alphaArm,pTArm,cosThStar,asym);
1662   //  }
1663 }
1664
1665 //_______________________________________________________________________________________
1666 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
1667                                               Float_t pt1,   Float_t pt2,
1668                                               Int_t ncell1,  Int_t ncell2,
1669                                               Double_t mass, Double_t pt,  Double_t asym,
1670                                               Double_t deta, Double_t dphi)
1671 {
1672   //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1673   //Adjusted for Pythia, need to see what to do for other generators.
1674   //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1675   // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1676   
1677   Int_t ancPDG    = 0;
1678   Int_t ancStatus = 0;
1679   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1680                                                               GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex);
1681   
1682   Int_t momindex  = -1;
1683   Int_t mompdg    = -1;
1684   Int_t momstatus = -1;
1685   
1686   Float_t prodR = -1;
1687   Int_t mcIndex = -1;
1688   
1689   if(ancLabel > -1)
1690   {
1691     AliDebug(1,Form("Common ancestor label %d, pdg %d, name %s, status %d",
1692                     ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus));
1693     
1694     if(ancPDG==22)
1695     {//gamma
1696       mcIndex = 0;
1697     }
1698     else if(TMath::Abs(ancPDG)==11)
1699     {//e
1700       mcIndex = 1;
1701     }
1702     else if(ancPDG==111)
1703     {//Pi0
1704       mcIndex = 2;
1705       if(fMultiCutAnaSim)
1706       {
1707         for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1708         {
1709           for(Int_t icell=0; icell<fNCellNCuts; icell++)
1710           {
1711             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1712             {
1713               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1714               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1715                  asym   <  fAsymCuts[iasym]                                   &&
1716                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1717               {
1718                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1719                 fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass);
1720                 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(),pt);
1721               }//pass the different cuts
1722             }// pid bit cut loop
1723           }// icell loop
1724         }// pt cut loop
1725       }//Multi cut ana sim
1726       else
1727       {
1728         fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(),mass);
1729         
1730         if(mass < 0.17 && mass > 0.1)
1731         {
1732           fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(),pt);
1733           
1734           //Int_t uniqueId = -1;
1735           if(GetReader()->ReadStack())
1736           {
1737             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1738             momindex  = ancestor->GetFirstMother();
1739             if(momindex < 0) return;
1740             TParticle* mother = GetMCStack()->Particle(momindex);
1741             mompdg    = TMath::Abs(mother->GetPdgCode());
1742             momstatus = mother->GetStatusCode();
1743             prodR = mother->R();
1744             //uniqueId = mother->GetUniqueID();
1745           }
1746           else
1747           {
1748             TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1749             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1750             momindex  = ancestor->GetMother();
1751             if(momindex < 0) return;
1752             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1753             mompdg    = TMath::Abs(mother->GetPdgCode());
1754             momstatus = mother->GetStatus();
1755             prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1756             //uniqueId = mother->GetUniqueID();
1757           }
1758           
1759           //            printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1760           //                   pt,prodR,mompdg,momstatus,uniqueId);
1761           
1762           fhMCPi0ProdVertex->Fill(pt,prodR);
1763           
1764           if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1765           else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1766           else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1767           else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1768           else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1769           else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1770           else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1771           else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1772           else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1773           else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1774           else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1775           
1776           
1777         }//pi0 mass region
1778       }
1779     }
1780     else if(ancPDG==221)
1781     {//Eta
1782       mcIndex = 3;
1783       if(fMultiCutAnaSim)
1784       {
1785         for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1786         {
1787           for(Int_t icell=0; icell<fNCellNCuts; icell++)
1788           {
1789             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1790             {
1791               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1792               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1793                  asym   <  fAsymCuts[iasym]                                   &&
1794                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1795               {
1796                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1797                 fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(),mass);
1798                 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(),pt);
1799               }//pass the different cuts
1800             }// pid bit cut loop
1801           }// icell loop
1802         }// pt cut loop
1803       } //Multi cut ana sim
1804       else
1805       {
1806         fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(),mass);
1807         if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.Pt(),pt);
1808         
1809         if(GetReader()->ReadStack())
1810         {
1811           TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1812           momindex  = ancestor->GetFirstMother();
1813           if(momindex < 0) return;
1814           TParticle* mother = GetMCStack()->Particle(momindex);
1815           mompdg    = TMath::Abs(mother->GetPdgCode());
1816           momstatus = mother->GetStatusCode();
1817           prodR = mother->R();
1818         }
1819         else
1820         {
1821           TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1822           AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1823           momindex  = ancestor->GetMother();
1824           if(momindex < 0) return;
1825           AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1826           mompdg    = TMath::Abs(mother->GetPdgCode());
1827           momstatus = mother->GetStatus();
1828           prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1829         }
1830         
1831         fhMCEtaProdVertex->Fill(pt,prodR);
1832         
1833         if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1834         else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1835         else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1836         else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1837         else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1838         else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1839         //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1840         
1841       }// eta mass region
1842     }
1843     else if(ancPDG==-2212){//AProton
1844       mcIndex = 4;
1845     }
1846     else if(ancPDG==-2112){//ANeutron
1847       mcIndex = 5;
1848     }
1849     else if(TMath::Abs(ancPDG)==13){//muons
1850       mcIndex = 6;
1851     }
1852     else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
1853     {
1854       if(ancStatus==1)
1855       {//Stable particles, converted? not decayed resonances
1856         mcIndex = 6;
1857       }
1858       else
1859       {//resonances and other decays, more hadron conversions?
1860         mcIndex = 7;
1861       }
1862     }
1863     else
1864     {//Partons, colliding protons, strings, intermediate corrections
1865       if(ancStatus==11 || ancStatus==12)
1866       {//String fragmentation
1867         mcIndex = 8;
1868       }
1869       else if (ancStatus==21){
1870         if(ancLabel < 2)
1871         {//Colliding protons
1872           mcIndex = 11;
1873         }//colliding protons
1874         else if(ancLabel < 6)
1875         {//partonic initial states interactions
1876           mcIndex = 9;
1877         }
1878         else if(ancLabel < 8)
1879         {//Final state partons radiations?
1880           mcIndex = 10;
1881         }
1882         // else {
1883         //   printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1884         //          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1885         // }
1886       }//status 21
1887       //else {
1888       //  printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1889       //         ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1890       // }
1891     }////Partons, colliding protons, strings, intermediate corrections
1892   }//ancLabel > -1
1893   else
1894   { //ancLabel <= -1
1895     //printf("Not related at all label = %d\n",ancLabel);
1896     AliDebug(1,"Common ancestor not found");
1897     
1898     mcIndex = 12;
1899   }
1900   
1901   if(mcIndex >= 0 && mcIndex < 13)
1902   {
1903     fhMCOrgMass[mcIndex]->Fill(pt,mass);
1904     fhMCOrgAsym[mcIndex]->Fill(pt,asym);
1905     fhMCOrgDeltaEta[mcIndex]->Fill(pt,deta);
1906     fhMCOrgDeltaPhi[mcIndex]->Fill(pt,dphi);
1907   }
1908   
1909 }
1910
1911 //__________________________________________
1912 void AliAnaPi0::MakeAnalysisFillHistograms()
1913 {
1914   //Process one event and extract photons from AOD branch
1915   // filled with AliAnaPhoton and fill histos with invariant mass
1916   
1917   //In case of simulated data, fill acceptance histograms
1918   if(IsDataMC())FillAcceptanceHistograms();
1919   
1920   //if (GetReader()->GetEventNumber()%10000 == 0)
1921   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1922   
1923   if(!GetInputAODBranch())
1924   {
1925     AliFatal(Form("No input aod photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
1926     return;
1927   }
1928   
1929   //Init some variables
1930   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
1931   
1932   AliDebug(1,Form("Photon entries %d", nPhot));
1933   
1934   //If less than photon 2 entries in the list, skip this event
1935   if ( nPhot < 2 )
1936   {
1937     AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentrality()));
1938     
1939     if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
1940     
1941     return ;
1942   }
1943   
1944   Int_t ncentr = GetNCentrBin();
1945   if(GetNCentrBin()==0) ncentr = 1;
1946   
1947   //Init variables
1948   Int_t module1         = -1;
1949   Int_t module2         = -1;
1950   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
1951   Int_t evtIndex1       = 0 ;
1952   Int_t currentEvtIndex = -1;
1953   Int_t curCentrBin     = GetEventCentralityBin();
1954   //Int_t curVzBin        = GetEventVzBin();
1955   //Int_t curRPBin        = GetEventRPBin();
1956   Int_t eventbin        = GetEventMixBin();
1957   
1958   if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
1959   {
1960     AliWarning(Form("Mix Bin not expected: cen bin %d, z bin %d, rp bin %d, total bin %d, Event Centrality %d, z vertex %2.3f, Reaction Plane %2.3f",
1961                     GetEventCentralityBin(),GetEventVzBin(), GetEventRPBin(),eventbin,GetEventCentrality(),vert[2],GetEventPlaneAngle()));
1962     return;
1963   }
1964   
1965   //Get shower shape information of clusters
1966   TObjArray *clusters = 0;
1967   if     (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
1968   else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ;
1969   
1970   //---------------------------------
1971   //First loop on photons/clusters
1972   //---------------------------------
1973   for(Int_t i1=0; i1<nPhot-1; i1++)
1974   {
1975     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1976     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
1977     
1978     // get the event index in the mixed buffer where the photon comes from
1979     // in case of mixing with analysis frame, not own mixing
1980     evtIndex1 = GetEventIndex(p1, vert) ;
1981     if ( evtIndex1 == -1 )
1982       return ;
1983     if ( evtIndex1 == -2 )
1984       continue ;
1985     
1986     // Only effective in case of mixed event frame
1987     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
1988     
1989     if (evtIndex1 != currentEvtIndex)
1990     {
1991       //Fill event bin info
1992       if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
1993       if(GetNCentrBin() > 1)
1994       {
1995         fhCentrality->Fill(curCentrBin);
1996         if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
1997       }
1998       currentEvtIndex = evtIndex1 ;
1999     }
2000     
2001     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2002     
2003     //Get the momentum of this cluster
2004     fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2005     
2006     //Get (Super)Module number of this cluster
2007     module1 = GetModuleNumber(p1);
2008     
2009     //------------------------------------------
2010     // Recover original cluster
2011     Int_t iclus1 = -1 ;
2012     AliVCluster * cluster1 = FindCluster(clusters,p1->GetCaloLabel(0),iclus1);
2013     if(!cluster1) AliWarning("Cluster1 not found!");
2014     
2015     //---------------------------------
2016     //Second loop on photons/clusters
2017     //---------------------------------
2018     for(Int_t i2=i1+1; i2<nPhot; i2++)
2019     {
2020       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2021       //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
2022       
2023       //In case of mixing frame, check we are not in the same event as the first cluster
2024       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2025       if ( evtIndex2 == -1 )
2026         return ;
2027       if ( evtIndex2 == -2 )
2028         continue ;
2029       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2030         continue ;
2031       
2032       //------------------------------------------
2033       // Recover original cluster
2034       Int_t iclus2 = -1;
2035       AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
2036       // start new loop from iclus1+1 to gain some time
2037       if(!cluster2) AliWarning("Cluster2 not found!");
2038       
2039       // Get the TOF,l0 and ncells from the clusters
2040       Float_t tof1  = -1;
2041       Float_t l01   = -1;
2042       Int_t ncell1  = 0;
2043       if(cluster1)
2044       {
2045         tof1   = cluster1->GetTOF()*1e9;
2046         l01    = cluster1->GetM02();
2047         ncell1 = cluster1->GetNCells();
2048         //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
2049       }
2050       //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2051       //            p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2052       
2053       Float_t tof2  = -1;
2054       Float_t l02   = -1;
2055       Int_t ncell2  = 0;
2056       if(cluster2)
2057       {
2058         tof2   = cluster2->GetTOF()*1e9;
2059         l02    = cluster2->GetM02();
2060         ncell2 = cluster2->GetNCells();
2061         //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
2062       }
2063       //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2064       //            p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2065       
2066       if(cluster1 && cluster2)
2067       {
2068         Double_t t12diff = tof1-tof2;
2069         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2070       }
2071       //------------------------------------------
2072       
2073       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2074       
2075       //Get the momentum of this cluster
2076       fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2077       
2078       //Get module number
2079       module2       = GetModuleNumber(p2);
2080       
2081       //---------------------------------
2082       // Get pair kinematics
2083       //---------------------------------
2084       Double_t m    = (fPhotonMom1 + fPhotonMom2).M() ;
2085       Double_t pt   = (fPhotonMom1 + fPhotonMom2).Pt();
2086       Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
2087       Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
2088       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2089       
2090       AliDebug(2,Form("E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a));
2091       
2092       //--------------------------------
2093       // Opening angle selection
2094       //--------------------------------
2095       //Check if opening angle is too large or too small compared to what is expected
2096       Double_t angle   = fPhotonMom1.Angle(fPhotonMom2.Vect());
2097       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
2098       {
2099         AliDebug(2,Form("Real pair angle %f not in E %f window",angle, (fPhotonMom1+fPhotonMom2).E()));
2100         continue;
2101       }
2102       
2103       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2104       {
2105         AliDebug(2,Form("Real pair cut %f < angle %f < cut %f",fAngleCut, angle, fAngleMaxCut));
2106         continue;
2107       }
2108       
2109       //-------------------------------------------------------------------------------------------------
2110       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2111       //-------------------------------------------------------------------------------------------------
2112       if(a < fAsymCuts[0] && fFillSMCombinations)
2113       {
2114         if(module1==module2 && module1 >=0 && module1<fNModules)
2115           fhReMod[module1]->Fill(pt,m) ;
2116         
2117         if(GetCalorimeter()==kEMCAL)
2118         {
2119           // Same sector
2120           Int_t j=0;
2121           for(Int_t i = 0; i < fNModules/2; i++)
2122           {
2123             j=2*i;
2124             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2125           }
2126           
2127           // Same side
2128           for(Int_t i = 0; i < fNModules-2; i++){
2129             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2130           }
2131         }//EMCAL
2132         else {//PHOS
2133           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2134           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2135           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2136         }//PHOS
2137       }
2138       
2139       //In case we want only pairs in same (super) module, check their origin.
2140       Bool_t ok = kTRUE;
2141       if(fSameSM && module1!=module2) ok=kFALSE;
2142       if(ok)
2143       {
2144         //Check if one of the clusters comes from a conversion
2145         if(fCheckConversion)
2146         {
2147           if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2148           else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2149         }
2150         
2151         // Fill shower shape cut histograms
2152         if(fFillSSCombinations)
2153         {
2154           if     ( l01 > 0.01 && l01 < 0.4  &&
2155                   l02 > 0.01 && l02 < 0.4 )               fhReSS[0]->Fill(pt,m); // Tight
2156           else if( l01 > 0.4  && l02 > 0.4 )               fhReSS[1]->Fill(pt,m); // Loose
2157           else if( l01 > 0.01 && l01 < 0.4  && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2158           else if( l02 > 0.01 && l02 < 0.4  && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2159         }
2160         
2161         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2162         for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2163         {
2164           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
2165           {
2166             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2167             {
2168               if(a < fAsymCuts[iasym])
2169               {
2170                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2171                 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d - max index %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym, curCentrBin*fNPIDBits*fNAsymCuts);
2172                 
2173                 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2174                 
2175                 fhRe1     [index]->Fill(pt,m);
2176                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2177                 if(fFillBadDistHisto)
2178                 {
2179                   if(p1->DistToBad()>0 && p2->DistToBad()>0)
2180                   {
2181                     fhRe2     [index]->Fill(pt,m) ;
2182                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2183                     if(p1->DistToBad()>1 && p2->DistToBad()>1)
2184                     {
2185                       fhRe3     [index]->Fill(pt,m) ;
2186                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2187                     }// bad 3
2188                   }// bad2
2189                 }// Fill bad dist histos
2190               }//assymetry cut
2191             }// asymmetry cut loop
2192           }// bad 1
2193         }// pid bit loop
2194         
2195         //Fill histograms with opening angle
2196         if(fFillAngleHisto)
2197         {
2198           fhRealOpeningAngle   ->Fill(pt,angle);
2199           fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2200         }
2201         
2202         //Fill histograms with pair assymmetry
2203         if(fFillAsymmetryHisto)
2204         {
2205           fhRePtAsym->Fill(pt,a);
2206           if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2207           if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2208         }
2209         
2210         //---------
2211         // MC data
2212         //---------
2213         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2214         if(IsDataMC())
2215         {
2216           if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2217              GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2218           {
2219             fhReMCFromConversion->Fill(pt,m);
2220           }
2221           else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2222                   !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2223           {
2224             fhReMCFromNotConversion->Fill(pt,m);
2225           }
2226           else
2227           {
2228             fhReMCFromMixConversion->Fill(pt,m);
2229           }
2230           
2231           if(fFillOriginHisto)
2232             FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2233         }
2234         
2235         //-----------------------
2236         //Multi cuts analysis
2237         //-----------------------
2238         if(fMultiCutAna)
2239         {
2240           //Histograms for different PID bits selection
2241           for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2242           {
2243             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
2244                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
2245             
2246             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2247           } // pid bit cut loop
2248           
2249           //Several pt,ncell and asymmetry cuts
2250           for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2251           {
2252             for(Int_t icell=0; icell<fNCellNCuts; icell++)
2253             {
2254               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2255               {
2256                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2257                 if(p1->E() >   fPtCuts[ipt]      && p2->E() > fPtCuts[ipt]        &&
2258                    a        <   fAsymCuts[iasym]                                  &&
2259                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell])
2260                 {
2261                   fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2262                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2263                   if(fFillSMCombinations && module1==module2)
2264                   {
2265                     fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2266                   }
2267                 }
2268               }// pid bit cut loop
2269             }// icell loop
2270           }// pt cut loop
2271           
2272           if(GetHistogramRanges()->GetHistoTrackMultiplicityBins())
2273           {
2274             for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2275             {
2276               if(a < fAsymCuts[iasym]) fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2277             }
2278           }
2279         }// multiple cuts analysis
2280         
2281       }// ok if same sm
2282       
2283     }// second same event particle
2284   }// first cluster
2285   
2286   //-------------------------------------------------------------
2287   // Mixing
2288   //-------------------------------------------------------------
2289   if(DoOwnMix())
2290   {
2291     //Recover events in with same characteristics as the current event
2292     
2293     //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2294     if(eventbin < 0) return ;
2295     
2296     TList * evMixList=fEventsList[eventbin] ;
2297     
2298     if(!evMixList)
2299     {
2300       AliWarning(Form("Mix event list not available, bin %d",eventbin));
2301       return;
2302     }
2303     
2304     Int_t nMixed = evMixList->GetSize() ;
2305     for(Int_t ii=0; ii<nMixed; ii++)
2306     {
2307       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2308       Int_t nPhot2=ev2->GetEntriesFast() ;
2309       Double_t m = -999;
2310       AliDebug(1,Form("Mixed event %d photon entries %d, centrality bin %d",ii, nPhot2, GetEventCentralityBin()));
2311       
2312       fhEventMixBin->Fill(eventbin) ;
2313       
2314       //---------------------------------
2315       //First loop on photons/clusters
2316       //---------------------------------
2317       for(Int_t i1=0; i1<nPhot; i1++)
2318       {
2319         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2320         
2321         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2322         
2323         //Get kinematics of cluster and (super) module of this cluster
2324         fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2325         module1 = GetModuleNumber(p1);
2326         
2327         //---------------------------------
2328         //First loop on photons/clusters
2329         //---------------------------------
2330         for(Int_t i2=0; i2<nPhot2; i2++)
2331         {
2332           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2333           
2334           //Get kinematics of second cluster and calculate those of the pair
2335           fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2336           m           = (fPhotonMom1+fPhotonMom2).M() ;
2337           Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2338           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2339           
2340           //Check if opening angle is too large or too small compared to what is expected
2341           Double_t angle   = fPhotonMom1.Angle(fPhotonMom2.Vect());
2342           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
2343           {
2344             AliDebug(2,Form("Mix pair angle %f not in E %f window",angle, (fPhotonMom1+fPhotonMom2).E()));
2345             continue;
2346           }
2347           
2348           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2349           {
2350             AliDebug(2,Form("Mix pair angle %f < cut %f",angle,fAngleCut));
2351             continue;
2352           }
2353           
2354           AliDebug(2,Form("Mixed Event: pT: fPhotonMom1 %2.2f, fPhotonMom2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %2.3f",p1->Pt(), p2->Pt(), pt,m,a));
2355           
2356           //In case we want only pairs in same (super) module, check their origin.
2357           module2 = GetModuleNumber(p2);
2358           
2359           //-------------------------------------------------------------------------------------------------
2360           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2361           //-------------------------------------------------------------------------------------------------
2362           if(a < fAsymCuts[0] && fFillSMCombinations){
2363             if(module1==module2 && module1 >=0 && module1<fNModules)
2364               fhMiMod[module1]->Fill(pt,m) ;
2365             
2366             if(GetCalorimeter()==kEMCAL)
2367             {
2368               // Same sector
2369               Int_t j=0;
2370               for(Int_t i = 0; i < fNModules/2; i++)
2371               {
2372                 j=2*i;
2373                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2374               }
2375               
2376               // Same side
2377               for(Int_t i = 0; i < fNModules-2; i++)
2378               {
2379                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2380               }
2381             }//EMCAL
2382             else
2383             {//PHOS
2384               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2385               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2386               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2387             }//PHOS
2388             
2389             
2390           }
2391           
2392           Bool_t ok = kTRUE;
2393           if(fSameSM && module1!=module2) ok=kFALSE;
2394           if(ok){
2395             
2396             //Check if one of the clusters comes from a conversion
2397             if(fCheckConversion)
2398             {
2399               if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2400               else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2401             }
2402             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2403             for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2404             {
2405               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2406               {
2407                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2408                 {
2409                   if(a < fAsymCuts[iasym])
2410                   {
2411                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2412                     
2413                     if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2414                     
2415                     fhMi1     [index]->Fill(pt,m) ;
2416                     
2417                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2418                     
2419                     if(fFillBadDistHisto)
2420                     {
2421                       if(p1->DistToBad()>0 && p2->DistToBad()>0)
2422                       {
2423                         fhMi2     [index]->Fill(pt,m) ;
2424                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2425                         if(p1->DistToBad()>1 && p2->DistToBad()>1)
2426                         {
2427                           fhMi3     [index]->Fill(pt,m) ;
2428                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2429                         }
2430                       }
2431                     }// Fill bad dist histo
2432                   }//Asymmetry cut
2433                 }// Asymmetry loop
2434               }//PID cut
2435             }//loop for histograms
2436             
2437             //-----------------------
2438             //Multi cuts analysis
2439             //-----------------------
2440             if(fMultiCutAna){
2441               //Several pt,ncell and asymmetry cuts
2442               
2443               for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2444               {
2445                 for(Int_t icell=0; icell<fNCellNCuts; icell++)
2446                 {
2447                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2448                   {
2449                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2450                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2451                        a        <   fAsymCuts[iasym]                                    //&&
2452                        //p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2453                        )
2454                     {
2455                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2456                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2457                     }
2458                   }// pid bit cut loop
2459                 }// icell loop
2460               }// pt cut loop
2461             } // Multi cut ana
2462             
2463             //Fill histograms with opening angle
2464             if(fFillAngleHisto)
2465             {
2466               fhMixedOpeningAngle   ->Fill(pt,angle);
2467               fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2468             }
2469             
2470           }//ok
2471         }// second cluster loop
2472       }//first cluster loop
2473     }//loop on mixed events
2474     
2475     //--------------------------------------------------------
2476     // Add the current event to the list of events for mixing
2477     //--------------------------------------------------------
2478     
2479     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2480     //Add current event to buffer and Remove redundant events
2481     if( currentEvent->GetEntriesFast() > 0 )
2482     {
2483       evMixList->AddFirst(currentEvent) ;
2484       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2485       if( evMixList->GetSize() >= GetNMaxEvMix() )
2486       {
2487         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2488         evMixList->RemoveLast() ;
2489         delete tmp ;
2490       }
2491     }
2492     else
2493     { //empty event
2494       delete currentEvent ;
2495       currentEvent=0 ;
2496     }
2497   }// DoOwnMix
2498   
2499   AliDebug(1,"End fill histograms");
2500 }
2501
2502 //________________________________________________________________________
2503 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2504 {
2505   // retieves the event index and checks the vertex
2506   //    in the mixed buffer returns -2 if vertex NOK
2507   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2508   
2509   Int_t evtIndex = -1 ;
2510   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2511   {
2512     if (GetMixedEvent())
2513     {
2514       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2515       GetVertex(vert,evtIndex);
2516       
2517       if(TMath::Abs(vert[2])> GetZvertexCut())
2518         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2519     }
2520     else
2521     {
2522       // Single event
2523       GetVertex(vert);
2524       
2525       if(TMath::Abs(vert[2])> GetZvertexCut())
2526         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2527       else
2528         evtIndex = 0 ;
2529     }
2530   }//No MC reader
2531   else
2532   {
2533     evtIndex = 0;
2534     vert[0] = 0. ; 
2535     vert[1] = 0. ; 
2536     vert[2] = 0. ; 
2537   }
2538   
2539   return evtIndex ; 
2540 }
2541