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