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