fix the selection of the meson objects when they are identified as pi0 or eta, before...
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPi0EbE.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 for the analysis of high pT pi0 event by event
18 // Pi0/Eta identified by one of the following:
19 //  -Invariant mass of 2 cluster in calorimeter
20 //  -Shower shape analysis in calorimeter
21 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
22 //
23 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
25
26
27 // --- ROOT system ---
28 #include <TList.h>
29 #include <TClonesArray.h>
30 #include <TObjString.h>
31
32 // --- Analysis system ---
33 #include "AliAnaPi0EbE.h"
34 #include "AliCaloTrackReader.h"
35 #include "AliIsolationCut.h"
36 #include "AliNeutralMesonSelection.h"
37 #include "AliCaloPID.h"
38 #include "AliMCAnalysisUtils.h"
39 #include "AliStack.h"
40 #include "AliFiducialCut.h"
41 #include "TParticle.h"
42 #include "AliVCluster.h"
43 #include "AliESDEvent.h"
44 #include "AliAODEvent.h"
45 #include "AliAODMCParticle.h"
46
47 ClassImp(AliAnaPi0EbE)
48
49 //____________________________
50 AliAnaPi0EbE::AliAnaPi0EbE() :
51 AliAnaCaloTrackCorrBaseClass(),     fAnaType(kIMCalo),                  fCalorimeter(""),
52 fMinDist(0.),fMinDist2(0.),         fMinDist3(0.),
53 fNLMCutMin(-1),                     fNLMCutMax(10),
54 fTimeCutMin(-10000),                fTimeCutMax(10000),
55 fRejectTrackMatch(kTRUE),
56 fFillPileUpHistograms(0),
57 fFillWeightHistograms(kFALSE),      fFillTMHisto(0),
58 fFillSelectClHisto(0),              fFillOnlySimpleSSHisto(1),          fFillEMCALBCHistograms(0),
59 fInputAODGammaConvName(""),
60 fCheckSplitDistToBad(0),
61 // Histograms
62 fhPt(0),                            fhE(0),
63 fhPtEta(0),                         fhPtPhi(0),                         fhEtaPhi(0),
64 fhEtaPhiEMCALBC0(0),                fhEtaPhiEMCALBC1(0),                fhEtaPhiEMCALBCN(0),
65 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
66 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
67 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
68 fhPtCentrality(),                   fhPtEventPlane(0),
69 fhPtReject(0),                      fhEReject(0),
70 fhPtEtaReject(0),                   fhPtPhiReject(0),                   fhEtaPhiReject(0),
71 fhMass(0),                          fhMassPt(0),                        fhMassSplitPt(0),
72 fhSelectedMass(),                   fhSelectedMassPt(0),                fhSelectedMassSplitPt(0),
73 fhMassNoOverlap(0),                 fhMassPtNoOverlap(0),               fhMassSplitPtNoOverlap(0),
74 fhSelectedMassNoOverlap(0),         fhSelectedMassPtNoOverlap(0),       fhSelectedMassSplitPtNoOverlap(0),
75 fhMCPi0PtRecoPtPrim(0),                       fhMCEtaPtRecoPtPrim(0),
76 fhMCPi0PtRecoPtPrimNoOverlap(0),              fhMCEtaPtRecoPtPrimNoOverlap(0),
77 fhMCPi0SplitPtRecoPtPrim(0),                  fhMCEtaSplitPtRecoPtPrim(0),
78 fhMCPi0SplitPtRecoPtPrimNoOverlap(0),         fhMCEtaSplitPtRecoPtPrimNoOverlap(0),
79 fhMCPi0SelectedPtRecoPtPrim(0),               fhMCEtaSelectedPtRecoPtPrim(0),
80 fhMCPi0SelectedPtRecoPtPrimNoOverlap(0),      fhMCEtaSelectedPtRecoPtPrimNoOverlap(0),
81 fhMCPi0SelectedSplitPtRecoPtPrim(0),          fhMCEtaSelectedSplitPtRecoPtPrim(0),
82 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap(0),
83 fhAsymmetry(0),                     fhSelectedAsymmetry(0),
84 fhSplitE(0),                        fhSplitPt(0),
85 fhSplitPtEta(0),                    fhSplitPtPhi(0),
86 fhNLocMaxSplitPt(0),
87 fhPtDecay(0),                       fhEDecay(0),
88 // Shower shape histos
89 fhPtDispersion(0),                  fhPtLambda0(0),                     fhPtLambda1(0),
90 fhPtLambda0NoTRD(0),                fhPtLambda0FracMaxCellCut(0),
91 fhPtFracMaxCell(0),                 fhPtFracMaxCellNoTRD(0),
92 fhPtNCells(0),                      fhPtTime(0),                        fhEPairDiffTime(0),
93 fhPtDispEta(0),                     fhPtDispPhi(0),
94 fhPtSumEta(0),                      fhPtSumPhi(0),                      fhPtSumEtaPhi(0),
95 fhPtDispEtaPhiDiff(0),              fhPtSphericity(0),
96
97 // MC histos
98 fhMCE(),                            fhMCPt(),
99 fhMCPtPhi(),                        fhMCPtEta(),
100 fhMCEReject(),                      fhMCPtReject(),
101 fhMCPtCentrality(),
102 fhMCPi0PtGenRecoFraction(0),        fhMCEtaPtGenRecoFraction(0),
103 fhMCPi0DecayPt(0),                  fhMCPi0DecayPtFraction(0),
104 fhMCEtaDecayPt(0),                  fhMCEtaDecayPtFraction(0),
105 fhMCOtherDecayPt(0),
106 fhMassPairMCPi0(0),                 fhMassPairMCEta(0),
107 fhAnglePairMCPi0(0),                fhAnglePairMCEta(0),
108 fhMCPi0PtOrigin(0x0),               fhMCEtaPtOrigin(0x0),
109 fhMCPi0ProdVertex(0),               fhMCEtaProdVertex(0),
110
111 // Weight studies
112 fhECellClusterRatio(0),             fhECellClusterLogRatio(0),
113 fhEMaxCellClusterRatio(0),          fhEMaxCellClusterLogRatio(0),
114 fhTrackMatchedDEta(0),              fhTrackMatchedDPhi(0),              fhTrackMatchedDEtaDPhi(0),
115 fhTrackMatchedDEtaPos(0),           fhTrackMatchedDPhiPos(0),           fhTrackMatchedDEtaDPhiPos(0),
116 fhTrackMatchedDEtaNeg(0),           fhTrackMatchedDPhiNeg(0),           fhTrackMatchedDEtaDPhiNeg(0),
117 fhTrackMatchedMCParticlePt(0),
118 fhTrackMatchedMCParticleDEta(0),    fhTrackMatchedMCParticleDPhi(0),
119 fhdEdx(0),                          fhEOverP(0),                        fhEOverPNoTRD(0),
120 // Number of local maxima in cluster
121 fhNLocMaxPt(0),                     fhNLocMaxPtReject(0),
122 // PileUp
123 fhTimePtNoCut(0),                   fhTimePtSPD(0),                     fhTimePtSPDMulti(0),
124 fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
125 fhTimeNPileUpVertContributors(0),
126 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
127 fhPtNPileUpSPDVtx(0),               fhPtNPileUpTrkVtx(0),
128 fhPtNPileUpSPDVtxTimeCut(0),        fhPtNPileUpTrkVtxTimeCut(0),
129 fhPtNPileUpSPDVtxTimeCut2(0),       fhPtNPileUpTrkVtxTimeCut2(0)
130 {
131   //default ctor
132   
133   for(Int_t i = 0; i < 6; i++)
134   {
135     fhMCE              [i] = 0;
136     fhMCPt             [i] = 0;
137     fhMCPtPhi          [i] = 0;
138     fhMCPtEta          [i] = 0;
139     fhMCPtCentrality   [i] = 0;
140     
141     fhMCSplitE         [i] = 0;
142     fhMCSplitPt        [i] = 0;
143     fhMCSplitPtPhi     [i] = 0;
144     fhMCSplitPtEta     [i] = 0;
145     
146     fhMCNLocMaxPt      [i] = 0;
147     fhMCNLocMaxSplitPt [i] = 0;
148     fhMCNLocMaxPtReject[i] = 0;
149     
150     fhMCPtLambda0       [i] = 0;
151     fhMCPtLambda0NoTRD  [i] = 0;
152     fhMCPtLambda0FracMaxCellCut[i]= 0;
153     fhMCPtFracMaxCell   [i] = 0;
154     fhMCPtLambda1       [i] = 0;
155     fhMCPtDispersion    [i] = 0;
156     
157     fhMCPtDispEta       [i] = 0;
158     fhMCPtDispPhi       [i] = 0;
159     fhMCPtSumEtaPhi     [i] = 0;
160     fhMCPtDispEtaPhiDiff[i] = 0;
161     fhMCPtSphericity    [i] = 0;
162     fhMCPtAsymmetry     [i] = 0;
163     
164     fhMCMassPt             [i]=0;
165     fhMCMassSplitPt        [i]=0;
166     fhMCSelectedMassPt     [i]=0;
167     fhMCSelectedMassSplitPt[i]=0;
168     
169     fhMCMassPtNoOverlap             [i]=0;
170     fhMCMassSplitPtNoOverlap        [i]=0;
171     fhMCSelectedMassPtNoOverlap     [i]=0;
172     fhMCSelectedMassSplitPtNoOverlap[i]=0;
173     
174     for(Int_t j = 0; j < 7; j++)
175     {
176       fhMCLambda0DispEta    [j][i] = 0;
177       fhMCLambda0DispPhi    [j][i] = 0;
178       fhMCDispEtaDispPhi    [j][i] = 0;
179       fhMCAsymmetryLambda0  [j][i] = 0;
180       fhMCAsymmetryDispEta  [j][i] = 0;
181       fhMCAsymmetryDispPhi  [j][i] = 0;
182     }
183   }
184   
185   for(Int_t j = 0; j < 7; j++)
186   {
187     fhLambda0DispEta    [j] = 0;
188     fhLambda0DispPhi    [j] = 0;
189     fhDispEtaDispPhi    [j] = 0;
190     fhAsymmetryLambda0  [j] = 0;
191     fhAsymmetryDispEta  [j] = 0;
192     fhAsymmetryDispPhi  [j] = 0;
193     
194     fhPtPileUp       [j] = 0;
195   }
196   
197   for(Int_t i = 0; i < 3; i++)
198   {
199     fhPtLambda0LocMax       [i] = 0;
200     fhPtLambda1LocMax       [i] = 0;
201     fhPtDispersionLocMax    [i] = 0;
202     fhPtDispEtaLocMax       [i] = 0;
203     fhPtDispPhiLocMax       [i] = 0;
204     fhPtSumEtaPhiLocMax     [i] = 0;
205     fhPtDispEtaPhiDiffLocMax[i] = 0;
206     fhPtSphericityLocMax    [i] = 0;
207     fhPtAsymmetryLocMax     [i] = 0;
208     fhMassPtLocMax          [i] = 0;
209     fhSelectedMassPtLocMax  [i] = 0;
210     for(Int_t ipart = 0; ipart<6; ipart++)
211     {
212       fhMCPtLambda0LocMax     [ipart][i] = 0;
213       fhMCSelectedMassPtLocMax[ipart][i] = 0;
214     }
215     
216     fhMCPi0PtRecoPtPrimLocMax             [i] = 0;
217     fhMCEtaPtRecoPtPrimLocMax             [i] = 0;
218     fhMCPi0SplitPtRecoPtPrimLocMax        [i] = 0;
219     fhMCEtaSplitPtRecoPtPrimLocMax        [i] = 0;
220
221     fhMCPi0SelectedPtRecoPtPrimLocMax     [i] = 0;
222     fhMCEtaSelectedPtRecoPtPrimLocMax     [i] = 0;
223     fhMCPi0SelectedSplitPtRecoPtPrimLocMax[i] = 0;
224     fhMCEtaSelectedSplitPtRecoPtPrimLocMax[i] = 0;
225
226   }
227   
228   //Weight studies
229   for(Int_t i =0; i < 14; i++){
230     fhLambda0ForW0[i] = 0;
231     //fhLambda1ForW0[i] = 0;
232     if(i<8)fhMassPairLocMax[i] = 0;
233   }
234   
235   for(Int_t i = 0; i < 11; i++)
236   {
237     fhEtaPhiTriggerEMCALBC       [i] = 0 ;
238     fhTimeTriggerEMCALBC         [i] = 0 ;
239     fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
240     
241     fhEtaPhiTriggerEMCALBCUM     [i] = 0 ;
242     fhTimeTriggerEMCALBCUM       [i] = 0 ;
243     
244   }
245   
246   for(Int_t iSM = 0; iSM < 22; iSM++)
247   {
248     fhNLocMaxPtSM[iSM] = 0;
249     for(Int_t inlm = 0; inlm < 3; inlm++)
250     {
251       fhSelectedMassPtLocMaxSM    [inlm][iSM] = 0;
252       fhSelectedLambda0PtLocMaxSM [inlm][iSM] = 0;
253     }
254   }
255   //Initialize parameters
256   InitParameters();
257   
258 }
259
260 //___________________________________________________________________________________
261 void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
262 {
263   // Fill some histograms to understand pile-up
264   if(!fFillPileUpHistograms) return;
265   
266   //printf("E %f, time %f\n",energy,time);
267   AliVEvent * event = GetReader()->GetInputEvent();
268   
269   fhTimePtNoCut->Fill(pt,time);
270   if(GetReader()->IsPileUpFromSPD())     
271   
272   if(GetReader()->IsPileUpFromSPD())             { fhPtPileUp[0]->Fill(pt); fhTimePtSPD     ->Fill(pt,time); }
273   if(GetReader()->IsPileUpFromEMCal())             fhPtPileUp[1]->Fill(pt);
274   if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPileUp[2]->Fill(pt);
275   if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPileUp[3]->Fill(pt);
276   if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPileUp[4]->Fill(pt);
277   if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPileUp[5]->Fill(pt);
278   if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
279   
280   if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
281   
282   // cells in cluster
283   
284   AliVCaloCells* cells = 0;
285   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
286   else                        cells = GetPHOSCells();
287
288   Float_t maxCellFraction = 0.;
289   Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
290   
291   Double_t tmax  = cells->GetCellTime(absIdMax);
292   GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
293   tmax*=1.e9;
294     
295   //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
296   if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
297   {
298     for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
299     {
300       Int_t absId  = calo->GetCellsAbsId()[ipos];
301       
302       if( absId == absIdMax ) continue ;
303       
304       Double_t timecell  = cells->GetCellTime(absId);
305       Float_t  amp       = cells->GetCellAmplitude(absId);
306       Int_t    bc        = GetReader()->GetInputEvent()->GetBunchCrossNumber();
307       GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
308       timecell*=1e9;
309       
310       Float_t diff = (tmax-timecell);
311             
312       if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
313       
314       if(GetReader()->IsPileUpFromSPD())
315       {
316         fhPtCellTimePileUp[0]->Fill(pt, timecell);
317         fhPtTimeDiffPileUp[0]->Fill(pt, diff);
318        }
319       
320       if(GetReader()->IsPileUpFromEMCal())
321       {
322         fhPtCellTimePileUp[1]->Fill(pt, timecell);
323         fhPtTimeDiffPileUp[1]->Fill(pt, diff);
324       }
325       
326       if(GetReader()->IsPileUpFromSPDOrEMCal())
327       {
328         fhPtCellTimePileUp[2]->Fill(pt, timecell);
329         fhPtTimeDiffPileUp[2]->Fill(pt, diff);
330       }
331       
332       if(GetReader()->IsPileUpFromSPDAndEMCal())
333       {
334         fhPtCellTimePileUp[3]->Fill(pt, timecell);
335         fhPtTimeDiffPileUp[3]->Fill(pt, diff);
336       }
337       
338       if(GetReader()->IsPileUpFromSPDAndNotEMCal())
339       {
340         fhPtCellTimePileUp[4]->Fill(pt, timecell);
341         fhPtTimeDiffPileUp[4]->Fill(pt, diff);
342       }
343       
344       if(GetReader()->IsPileUpFromEMCalAndNotSPD())
345       {
346         fhPtCellTimePileUp[5]->Fill(pt, timecell);
347         fhPtTimeDiffPileUp[5]->Fill(pt, diff);
348       }
349       
350       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
351       {
352         fhPtCellTimePileUp[6]->Fill(pt, timecell);
353         fhPtTimeDiffPileUp[6]->Fill(pt, diff);
354       }
355     }//loop
356   }
357
358   if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
359   
360   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
361   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
362   
363   // N pile up vertices
364   Int_t nVtxSPD = -1;
365   Int_t nVtxTrk = -1;
366   
367   if      (esdEv)
368   {
369     nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
370     nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
371     
372   }//ESD
373   else if (aodEv)
374   {
375     nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
376     nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
377   }//AOD
378   
379   fhTimeNPileUpVertSPD  ->Fill(time,nVtxSPD);
380   fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
381   
382         fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
383         fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
384         
385         if(TMath::Abs(time) < 25)
386         {
387                 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
388                 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
389   }
390   
391   if(time < 75 && time > -25)
392   {
393     fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
394     fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
395   }
396   
397   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
398   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
399   
400   Int_t ncont = -1;
401   Float_t z1 = -1, z2 = -1;
402   Float_t diamZ = -1;
403   for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
404   {
405     if      (esdEv)
406     {
407       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
408       ncont=pv->GetNContributors();
409       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
410       z2 = pv->GetZ();
411       diamZ = esdEv->GetDiamondZ();
412     }//ESD
413     else if (aodEv)
414     {
415       AliAODVertex *pv=aodEv->GetVertex(iVert);
416       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
417       ncont=pv->GetNContributors();
418       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
419       z2=pv->GetZ();
420       diamZ = aodEv->GetDiamondZ();
421     }// AOD
422     
423     Double_t distZ  = TMath::Abs(z2-z1);
424     diamZ  = TMath::Abs(z2-diamZ);
425     
426     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
427     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
428     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
429     
430   }// vertex loop
431 }
432
433
434 //______________________________________________________________________________________________
435 void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
436 {
437   // Fill histograms that do not pass the identification (SS case only)
438   
439   Float_t ener  = mom.E();
440   Float_t pt    = mom.Pt();
441   Float_t phi   = mom.Phi();
442   if(phi < 0) phi+=TMath::TwoPi();
443   Float_t eta = mom.Eta();
444   
445   fhPtReject     ->Fill(pt);
446   fhEReject      ->Fill(ener);
447   
448   fhPtEtaReject  ->Fill(ener,eta);
449   fhPtPhiReject  ->Fill(ener,phi);
450   fhEtaPhiReject ->Fill(eta,phi);
451   
452   fhNLocMaxPtReject->Fill(pt,nMaxima);
453
454   if(IsDataMC())
455   {
456     Int_t mcIndex = GetMCIndex(mctag);
457     fhMCEReject  [mcIndex] ->Fill(ener);
458     fhMCPtReject [mcIndex] ->Fill(pt);
459     fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
460   }
461 }
462
463 //___________________________________________________________________________________
464 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt, Int_t nMaxima,
465                                                  Int_t tag, Float_t asy)
466 {
467   // Fill shower shape, timing and other histograms for selected clusters from decay
468   
469   Float_t ener = cluster->E();
470   Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
471   Float_t l0   = cluster->GetM02();
472   Float_t l1   = cluster->GetM20();
473   Int_t   nSM  = GetModuleNumber(cluster);
474   
475   Int_t ptbin = -1;
476   if      (pt < 2 ) ptbin = 0;
477   else if (pt < 4 ) ptbin = 1;
478   else if (pt < 6 ) ptbin = 2;
479   else if (pt < 10) ptbin = 3;
480   else if (pt < 15) ptbin = 4;
481   else if (pt < 20) ptbin = 5;
482   else              ptbin = 6;
483   
484   Int_t indexMax = -1;
485   if     (nMaxima==1) indexMax = 0 ;
486   else if(nMaxima==2) indexMax = 1 ;
487   else                indexMax = 2 ;
488   
489   
490   AliVCaloCells * cell = 0x0;
491   if(fCalorimeter == "PHOS")
492     cell = GetPHOSCells();
493   else
494     cell = GetEMCALCells();
495   
496   Float_t maxCellFraction = 0;
497   GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
498   fhPtFracMaxCell->Fill(pt,maxCellFraction);
499   
500   FillWeightHistograms(cluster);
501   
502   fhPtDispersion->Fill(pt, disp);
503   fhPtLambda0   ->Fill(pt, l0  );
504   fhPtLambda1   ->Fill(pt, l1  );
505   
506   Float_t ll0  = 0., ll1  = 0.;
507   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
508   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
509   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
510   {
511     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
512                                                                                  ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
513     
514     fhPtDispEta       -> Fill(pt,dEta);
515     fhPtDispPhi       -> Fill(pt,dPhi);
516     fhPtSumEta        -> Fill(pt,sEta);
517     fhPtSumPhi        -> Fill(pt,sPhi);
518     fhPtSumEtaPhi     -> Fill(pt,sEtaPhi);
519     fhPtDispEtaPhiDiff-> Fill(pt,dPhi-dEta);
520     if(dEta+dPhi>0)fhPtSphericity-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
521     
522     fhDispEtaDispPhi[ptbin]->Fill(dEta,dPhi);
523     fhLambda0DispEta[ptbin]->Fill(l0  ,dEta);
524     fhLambda0DispPhi[ptbin]->Fill(l0  ,dPhi);
525     
526     if (fAnaType==kSSCalo)
527     {
528       // Asymmetry histograms
529       fhAsymmetryLambda0[ptbin]->Fill(l0 ,asy);
530       fhAsymmetryDispEta[ptbin]->Fill(dEta,asy);
531       fhAsymmetryDispPhi[ptbin]->Fill(dPhi,asy);
532     }
533   }
534   
535   fhNLocMaxPt->Fill(pt,nMaxima);
536   
537   if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
538     fhNLocMaxPtSM[nSM]->Fill(pt,nMaxima);
539   
540   fhPtLambda0LocMax   [indexMax]->Fill(pt,l0);
541   fhPtLambda1LocMax   [indexMax]->Fill(pt,l1);
542   fhPtDispersionLocMax[indexMax]->Fill(pt,disp);
543
544   if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
545   {
546     fhPtDispEtaLocMax       [indexMax]-> Fill(pt,dEta);
547     fhPtDispPhiLocMax       [indexMax]-> Fill(pt,dPhi);
548     fhPtSumEtaPhiLocMax     [indexMax]-> Fill(pt,sEtaPhi);
549     fhPtDispEtaPhiDiffLocMax[indexMax]-> Fill(pt,dPhi-dEta);
550     if(dEta+dPhi>0)       fhPtSphericityLocMax[indexMax]->Fill(pt,(dPhi-dEta)/(dEta+dPhi));
551     if(fAnaType==kSSCalo) fhPtAsymmetryLocMax [indexMax]->Fill(pt  ,asy);
552     
553   }
554   
555   if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >= 0 &&
556      GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
557   {
558     fhPtLambda0NoTRD    ->Fill(pt, l0  );
559     fhPtFracMaxCellNoTRD->Fill(pt,maxCellFraction);
560   }
561   
562   if(maxCellFraction < 0.5)
563     fhPtLambda0FracMaxCellCut->Fill(pt, l0  );
564   
565   fhPtTime  ->Fill(pt, cluster->GetTOF()*1.e9);
566   fhPtNCells->Fill(pt, cluster->GetNCells());
567   
568   // Fill Track matching control histograms
569   if(fFillTMHisto)
570   {
571     Float_t dZ  = cluster->GetTrackDz();
572     Float_t dR  = cluster->GetTrackDx();
573     
574     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
575     {
576       dR = 2000., dZ = 2000.;
577       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
578     }
579     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
580     
581     AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
582     
583     Bool_t positive = kFALSE;
584     if(track) positive = (track->Charge()>0);
585     
586     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
587     {
588       fhTrackMatchedDEta->Fill(pt,dZ);
589       fhTrackMatchedDPhi->Fill(pt,dR);
590       if(ener > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
591       
592       if(track)
593       {
594         if(positive)
595         {
596           fhTrackMatchedDEtaPos->Fill(pt,dZ);
597           fhTrackMatchedDPhiPos->Fill(pt,dR);
598           if(ener > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
599         }
600         else
601         {
602           fhTrackMatchedDEtaNeg->Fill(pt,dZ);
603           fhTrackMatchedDPhiNeg->Fill(pt,dR);
604           if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
605         }
606     }
607     }
608     // Check dEdx and E/p of matched clusters
609     
610     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
611     {      
612       if(track)
613       {
614         Float_t dEdx = track->GetTPCsignal();
615         fhdEdx->Fill(pt, dEdx);
616         
617         Float_t eOverp = cluster->E()/track->P();
618         fhEOverP->Fill(pt,  eOverp);
619         
620         // Change nSM for year > 2011 (< 4 in 2012-13, none after)
621         if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >= 0 &&
622            GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
623           fhEOverPNoTRD->Fill(pt,  eOverp);
624         
625       }
626       //else
627       //  printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
628       
629       if(IsDataMC())
630       {
631         Float_t mctag = -1;
632         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
633         {
634           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
635                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  2.5 ;
636           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  0.5 ;
637           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  1.5 ;
638           else                                                                                 mctag =  3.5 ;
639           
640         }
641         else
642         {
643           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
644                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  6.5 ;
645           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  4.5 ;
646           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  5.5 ;
647           else                                                                                 mctag =  7.5 ;
648         }
649         
650         fhTrackMatchedMCParticlePt   ->Fill(pt, mctag);
651         fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
652         fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
653         
654       }  // MC
655     }
656   }// Track matching histograms
657   
658   if(IsDataMC())
659   {
660     Int_t mcIndex = GetMCIndex(tag);
661     
662     fhMCPtLambda0[mcIndex]    ->Fill(pt, l0);
663     fhMCPtLambda1[mcIndex]    ->Fill(pt, l1);
664     fhMCPtDispersion[mcIndex] ->Fill(pt, disp);
665     fhMCPtFracMaxCell[mcIndex]->Fill(pt,maxCellFraction);
666     
667     fhMCPtLambda0LocMax     [mcIndex][indexMax]->Fill(pt,l0);
668
669     if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
670        GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
671       fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0  );
672     
673     if(maxCellFraction < 0.5)
674       fhMCPtLambda0FracMaxCellCut[mcIndex]->Fill(pt, l0  );
675     
676     if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
677     {
678       fhMCPtDispEta        [mcIndex]-> Fill(pt,dEta);
679       fhMCPtDispPhi        [mcIndex]-> Fill(pt,dPhi);
680       fhMCPtSumEtaPhi      [mcIndex]-> Fill(pt,sEtaPhi);
681       fhMCPtDispEtaPhiDiff [mcIndex]-> Fill(pt,dPhi-dEta);
682       if(dEta+dPhi>0)fhMCPtSphericity[mcIndex]-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
683       
684       if (fAnaType==kSSCalo)
685       {
686         fhMCAsymmetryLambda0[ptbin][mcIndex]->Fill(l0 ,asy);
687         fhMCAsymmetryDispEta[ptbin][mcIndex]->Fill(dEta,asy);
688         fhMCAsymmetryDispPhi[ptbin][mcIndex]->Fill(dPhi,asy);
689       }
690       
691       fhMCDispEtaDispPhi[ptbin][mcIndex]->Fill(dEta,dPhi);
692       fhMCLambda0DispEta[ptbin][mcIndex]->Fill(l0  ,dEta);
693       fhMCLambda0DispPhi[ptbin][mcIndex]->Fill(l0  ,dPhi);
694       
695     }
696     
697   }//MC
698   
699 }
700
701 //________________________________________________________
702 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
703 {
704   // Calculate weights and fill histograms
705   
706   if(!fFillWeightHistograms || GetMixedEvent()) return;
707   
708   AliVCaloCells* cells = 0;
709   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
710   else                        cells = GetPHOSCells();
711   
712   // First recalculate energy in case non linearity was applied
713   Float_t  energy = 0;
714   Float_t  ampMax = 0;
715   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
716   {
717     
718     Int_t id       = clus->GetCellsAbsId()[ipos];
719     
720     //Recalibrate cell energy if needed
721     Float_t amp = cells->GetCellAmplitude(id);
722     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
723     
724     energy    += amp;
725     
726     if(amp> ampMax)
727       ampMax = amp;
728     
729   } // energy loop
730   
731   if(energy <=0 )
732   {
733     printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
734     return;
735   }
736   
737   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
738   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
739   
740   //Get the ratio and log ratio to all cells in cluster
741   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
742   {
743     Int_t id       = clus->GetCellsAbsId()[ipos];
744     
745     //Recalibrate cell energy if needed
746     Float_t amp = cells->GetCellAmplitude(id);
747     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
748     
749     fhECellClusterRatio   ->Fill(energy,amp/energy);
750     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
751   }
752   
753   //Recalculate shower shape for different W0
754   if(fCalorimeter=="EMCAL"){
755     
756     Float_t l0org = clus->GetM02();
757     Float_t l1org = clus->GetM20();
758     Float_t dorg  = clus->GetDispersion();
759     
760     for(Int_t iw = 0; iw < 14; iw++)
761     {
762       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
763       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
764       
765       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
766       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
767       
768     } // w0 loop
769     
770     // Set the original values back
771     clus->SetM02(l0org);
772     clus->SetM20(l1org);
773     clus->SetDispersion(dorg);
774     
775   }// EMCAL
776 }
777
778 //__________________________________________
779 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
780 {
781         //Save parameters used for analysis
782   TString parList ; //this will be list of parameters used for this analysis.
783   const Int_t buffersize = 255;
784   char onePar[buffersize] ;
785   
786   snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
787   parList+=onePar ;
788   snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
789   parList+=onePar ;
790   
791   if(fAnaType == kSSCalo)
792   {
793     snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
794     parList+=onePar ;
795     snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
796     parList+=onePar ;
797     snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
798     parList+=onePar ;
799     snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
800     parList+=onePar ;
801   }
802   
803   //Get parameters set in base class.
804   parList += GetBaseParametersList() ;
805   
806   //Get parameters set in PID class.
807   if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
808   
809   return new TObjString(parList) ;
810 }
811
812 //_____________________________________________
813 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
814 {
815   // Create histograms to be saved in output file and
816   // store them in outputContainer
817   TList * outputContainer = new TList() ;
818   outputContainer->SetName("Pi0EbEHistos") ;
819   
820   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
821   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
822   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
823   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
824   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
825   Int_t tbins    = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin   = GetHistogramRanges()->GetHistoTimeMin();
826   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax   = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin   = GetHistogramRanges()->GetHistoNClusterCellMin();
827   
828   Int_t   nmassbins   = GetHistogramRanges()->GetHistoMassBins();
829   Float_t massmin     = GetHistogramRanges()->GetHistoMassMin();
830   Float_t massmax     = GetHistogramRanges()->GetHistoMassMax();
831   
832   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
833   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
834   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
835   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
836   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
837   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
838   
839   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
840   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
841   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
842   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
843   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
844   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
845   
846   Int_t   ntimptbins  = GetHistogramRanges()->GetHistoTimeBins();
847   Float_t timemax     = GetHistogramRanges()->GetHistoTimeMax();
848   Float_t timemin     = GetHistogramRanges()->GetHistoTimeMin();
849   
850   TString nlm[]   = {"1 Local Maxima","2 Local Maxima", "NLM > 2"};
851   TString ptype[] = {"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
852   TString pname[] = {"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
853   Int_t   bin[]   = {0,2,4,6,10,15,20,100}; // energy bins
854   
855   fhPt  = new TH1F("hPt","Number of identified  #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
856   fhPt->SetYTitle("#it{N}");
857   fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
858   outputContainer->Add(fhPt) ;
859   
860   fhE  = new TH1F("hE","Number of identified  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
861   fhE->SetYTitle("#it{N}");
862   fhE->SetXTitle("#it{E} (GeV)");
863   outputContainer->Add(fhE) ;
864   
865   fhPtPhi  = new TH2F
866   ("hPtPhi","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
867   fhPtPhi->SetYTitle("#phi (rad)");
868   fhPtPhi->SetXTitle("#it{E} (GeV)");
869   outputContainer->Add(fhPtPhi) ;
870   
871   fhPtEta  = new TH2F
872   ("hPtEta","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
873   fhPtEta->SetYTitle("#eta");
874   fhPtEta->SetXTitle("#it{E} (GeV)");
875   outputContainer->Add(fhPtEta) ;
876   
877   fhEtaPhi  = new TH2F
878   ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
879   fhEtaPhi->SetYTitle("#phi (rad)");
880   fhEtaPhi->SetXTitle("#eta");
881   outputContainer->Add(fhEtaPhi) ;
882   
883   if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
884   {
885     fhEtaPhiEMCALBC0  = new TH2F
886     ("hEtaPhiEMCALBC0","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
887     fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
888     fhEtaPhiEMCALBC0->SetXTitle("#eta");
889     outputContainer->Add(fhEtaPhiEMCALBC0) ;
890     
891     fhEtaPhiEMCALBC1  = new TH2F
892     ("hEtaPhiEMCALBC1","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with 25 < |#it{t}| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
893     fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
894     fhEtaPhiEMCALBC1->SetXTitle("#eta");
895     outputContainer->Add(fhEtaPhiEMCALBC1) ;
896     
897     fhEtaPhiEMCALBCN  = new TH2F
898     ("hEtaPhiEMCALBCN","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
899     fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
900     fhEtaPhiEMCALBCN->SetXTitle("#eta");
901     outputContainer->Add(fhEtaPhiEMCALBCN) ;
902     
903     for(Int_t i = 0; i < 11; i++)
904     {
905       fhEtaPhiTriggerEMCALBC[i] = new TH2F
906       (Form("hEtaPhiTriggerEMCALBC%d",i-5),
907        Form("meson #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
908        netabins,etamin,etamax,nphibins,phimin,phimax);
909       fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
910       fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
911       outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
912       
913       fhTimeTriggerEMCALBC[i] = new TH2F
914       (Form("hTimeTriggerEMCALBC%d",i-5),
915        Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
916        nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
917       fhTimeTriggerEMCALBC[i]->SetXTitle("#it{E} (GeV)");
918       fhTimeTriggerEMCALBC[i]->SetYTitle("#it{t} (ns)");
919       outputContainer->Add(fhTimeTriggerEMCALBC[i]);
920       
921       fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
922       (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
923        Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
924        nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
925       fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("#it{E} (GeV)");
926       fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("#it{t} (ns)");
927       outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
928       
929       fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
930       (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
931        Form("meson #it{E} > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
932        netabins,etamin,etamax,nphibins,phimin,phimax);
933       fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
934       fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
935       outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
936       
937       fhTimeTriggerEMCALBCUM[i] = new TH2F
938       (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
939        Form("meson #it{t} vs #it{E}, unmatched trigger EMCAL-BC=%d",i-5),
940        nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
941       fhTimeTriggerEMCALBCUM[i]->SetXTitle("#it{E} (GeV)");
942       fhTimeTriggerEMCALBCUM[i]->SetYTitle("#it{t} (ns)");
943       outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
944       
945     }
946     
947     fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
948                                                       "cluster #it{t} vs #it{E} of clusters, no match, rematch open time",
949                                                       nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
950     fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("#it{E} (GeV)");
951     fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("#it{t} (ns)");
952     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
953     
954     
955     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
956                                                         "cluster #it{t} vs #it{E} of clusters, no match, rematch with neigbour parches",
957                                                         nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
958     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("#it{E} (GeV)");
959     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("#it{t} (ns)");
960     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
961     
962     fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
963                                                   "cluster #it{t} vs #it{E} of clusters, no match, rematch open time and neigbour",
964                                                   nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
965     fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("#it{E} (GeV)");
966     fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("#it{t} (ns)");
967     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
968     
969   }
970   
971   fhPtCentrality  = new TH2F("hPtCentrality","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
972   fhPtCentrality->SetYTitle("centrality");
973   fhPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
974   outputContainer->Add(fhPtCentrality) ;
975   
976   fhPtEventPlane  = new TH2F("hPtEventPlane","event plane angle vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
977   fhPtEventPlane->SetYTitle("Event plane angle (rad)");
978   fhPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
979   outputContainer->Add(fhPtEventPlane) ;
980   
981   if(fAnaType == kSSCalo)
982   {
983     fhPtReject  = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
984     fhPtReject->SetYTitle("#it{N}");
985     fhPtReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
986     outputContainer->Add(fhPtReject) ;
987     
988     fhEReject  = new TH1F("hEReject","Number of rejected as  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
989     fhEReject->SetYTitle("#it{N}");
990     fhEReject->SetXTitle("#it{E} (GeV)");
991     outputContainer->Add(fhEReject) ;
992     
993     fhPtPhiReject  = new TH2F
994     ("hPtPhiReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
995     fhPtPhiReject->SetYTitle("#phi (rad)");
996     fhPtPhiReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
997     outputContainer->Add(fhPtPhiReject) ;
998     
999     fhPtEtaReject  = new TH2F
1000     ("hPtEtaReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1001     fhPtEtaReject->SetYTitle("#eta");
1002     fhPtEtaReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1003     outputContainer->Add(fhPtEtaReject) ;
1004     
1005     fhEtaPhiReject  = new TH2F
1006     ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1007     fhEtaPhiReject->SetYTitle("#phi (rad)");
1008     fhEtaPhiReject->SetXTitle("#eta");
1009     outputContainer->Add(fhEtaPhiReject) ;
1010   }
1011   
1012   fhMass  = new TH2F
1013   ("hMass","all pairs #it{M}: #it{E} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1014   fhMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1015   fhMass->SetXTitle("#it{E} (GeV)");
1016   outputContainer->Add(fhMass) ;
1017   
1018   fhSelectedMass  = new TH2F
1019   ("hSelectedMass","Selected #pi^{0} (#eta) pairs #it{M}: E vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1020   fhSelectedMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1021   fhSelectedMass->SetXTitle("#it{E} (GeV)");
1022   outputContainer->Add(fhSelectedMass) ;
1023   
1024   if(fAnaType == kSSCalo)
1025   {
1026     
1027     fhMassPt  = new TH2F
1028     ("hMassPt","all pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1029     fhMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1030     fhMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1031     outputContainer->Add(fhMassPt) ;
1032     
1033     fhSelectedMassPt  = new TH2F
1034     ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1035     fhSelectedMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1036     fhSelectedMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1037     outputContainer->Add(fhSelectedMassPt) ;
1038     
1039     for(Int_t inlm = 0; inlm < 3; inlm++)
1040     {
1041       fhMassPtLocMax[inlm]  = new TH2F
1042       (Form("hMassPtNLocMax%d",inlm+1),Form("all pairs #it{M}: #it{p}_{T} vs #it{M} and NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1043       fhMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1044       fhMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1045       outputContainer->Add(fhMassPtLocMax[inlm]) ;
1046       
1047       fhSelectedMassPtLocMax[inlm]  = new TH2F
1048       (Form("hSelectedMassPtLocMax%d",inlm+1),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1049       fhSelectedMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1050       fhSelectedMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1051       outputContainer->Add(fhSelectedMassPtLocMax[inlm]) ;
1052       
1053       for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1054       {
1055         fhSelectedMassPtLocMaxSM[inlm][iSM]  = new TH2F
1056         (Form("hSelectedMassPtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1057         fhSelectedMassPtLocMaxSM[inlm][iSM]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1058         fhSelectedMassPtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1059         outputContainer->Add(fhSelectedMassPtLocMaxSM[inlm][iSM]) ;
1060
1061         fhSelectedLambda0PtLocMaxSM[inlm][iSM]  = new TH2F
1062         (Form("hSelectedLambda0PtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #lambda_{0}^{2}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1063         fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetYTitle("#lambda_{0}^{2}");
1064         fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1065         outputContainer->Add(fhSelectedLambda0PtLocMaxSM[inlm][iSM]) ;
1066       }
1067       
1068       if(IsDataMC())
1069       {
1070         for(Int_t ipart = 0; ipart < 6; ipart++)
1071         {
1072           fhMCSelectedMassPtLocMax[ipart][inlm]  = new TH2F
1073           (Form("hSelectedMassPtLocMax%d_MC%s",inlm+1,pname[ipart].Data()),
1074            Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s, %s",nlm[inlm].Data(),pname[ipart].Data()),
1075            nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1076           fhMCSelectedMassPtLocMax[ipart][inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1077           fhMCSelectedMassPtLocMax[ipart][inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1078           outputContainer->Add(fhMCSelectedMassPtLocMax[ipart][inlm]) ;
1079         }
1080       }
1081     }
1082     
1083     if(IsDataMC())
1084     {
1085       fhMassNoOverlap  = new TH2F
1086       ("hMassNoOverlap","all pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1087       fhMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1088       fhMassNoOverlap->SetXTitle("#it{E} (GeV)");
1089       outputContainer->Add(fhMassNoOverlap) ;
1090       
1091       fhSelectedMassNoOverlap  = new TH2F
1092       ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1093       fhSelectedMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1094       fhSelectedMassNoOverlap->SetXTitle("#it{E} (GeV)");
1095       outputContainer->Add(fhSelectedMassNoOverlap) ;
1096       
1097       fhMassPtNoOverlap  = new TH2F
1098       ("hMassPtNoOverlap","all pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1099       fhMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1100       fhMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1101       outputContainer->Add(fhMassPtNoOverlap) ;
1102       
1103       fhSelectedMassPtNoOverlap  = new TH2F
1104       ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1105       fhSelectedMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1106       fhSelectedMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1107       outputContainer->Add(fhSelectedMassPtNoOverlap) ;
1108     }
1109   }
1110   
1111   if(fAnaType != kSSCalo)
1112   {
1113     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1114     fhPtDecay->SetYTitle("#it{N}");
1115     fhPtDecay->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1116     outputContainer->Add(fhPtDecay) ;
1117     
1118     fhEDecay  = new TH1F("hEDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1119     fhEDecay->SetYTitle("#it{N}");
1120     fhEDecay->SetXTitle("#it{E} (GeV)");
1121     outputContainer->Add(fhEDecay) ;
1122   }
1123   
1124   ////////
1125   
1126   if( fFillSelectClHisto )
1127   {
1128     fhPtDispersion  = new TH2F
1129     ("hPtDispersion","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1130     fhPtDispersion->SetYTitle("D^{2}");
1131     fhPtDispersion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1132     outputContainer->Add(fhPtDispersion) ;
1133     
1134     fhPtLambda0  = new TH2F
1135     ("hPtLambda0","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1136     fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
1137     fhPtLambda0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1138     outputContainer->Add(fhPtLambda0) ;
1139     
1140     fhPtLambda1  = new TH2F
1141     ("hPtLambda1","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1142     fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
1143     fhPtLambda1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1144     outputContainer->Add(fhPtLambda1) ;
1145     
1146     fhPtLambda0FracMaxCellCut  = new TH2F
1147     ("hPtLambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1148     fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1149     fhPtLambda0FracMaxCellCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1150     outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
1151     
1152     fhPtFracMaxCell  = new TH2F
1153     ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1154     fhPtFracMaxCell->SetYTitle("Fraction");
1155     fhPtFracMaxCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1156     outputContainer->Add(fhPtFracMaxCell) ;
1157     
1158     if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >=0 )
1159     {
1160       fhPtLambda0NoTRD  = new TH2F
1161       ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1162       fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1163       fhPtLambda0NoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1164       outputContainer->Add(fhPtLambda0NoTRD) ;
1165       
1166       fhPtFracMaxCellNoTRD  = new TH2F
1167       ("hPtFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1168       fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
1169       fhPtFracMaxCellNoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1170       outputContainer->Add(fhPtFracMaxCellNoTRD) ;
1171       
1172       if(!fFillOnlySimpleSSHisto)
1173       {
1174         fhPtDispEta  = new TH2F ("hPtDispEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1175         fhPtDispEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1176         fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
1177         outputContainer->Add(fhPtDispEta);
1178         
1179         fhPtDispPhi  = new TH2F ("hPtDispPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1180         fhPtDispPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1181         fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
1182         outputContainer->Add(fhPtDispPhi);
1183         
1184         fhPtSumEta  = new TH2F ("hPtSumEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs #it{p}_{T}",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1185         fhPtSumEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1186         fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
1187         outputContainer->Add(fhPtSumEta);
1188         
1189         fhPtSumPhi  = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs #it{p}_{T}",
1190                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1191         fhPtSumPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1192         fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
1193         outputContainer->Add(fhPtSumPhi);
1194         
1195         fhPtSumEtaPhi  = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",
1196                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1197         fhPtSumEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1198         fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
1199         outputContainer->Add(fhPtSumEtaPhi);
1200         
1201         fhPtDispEtaPhiDiff  = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",
1202                                        nptbins,ptmin,ptmax,200, -10,10);
1203         fhPtDispEtaPhiDiff->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1204         fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1205         outputContainer->Add(fhPtDispEtaPhiDiff);
1206         
1207         fhPtSphericity  = new TH2F ("hPtSphericity","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs #it{p}_{T} (GeV/#it{c})",
1208                                    nptbins,ptmin,ptmax, 200, -1,1);
1209         fhPtSphericity->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1210         fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1211         outputContainer->Add(fhPtSphericity);
1212         
1213         for(Int_t i = 0; i < 7; i++)
1214         {
1215           fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1216                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1217           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1218           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1219           outputContainer->Add(fhDispEtaDispPhi[i]);
1220           
1221           fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1222                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1223           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1224           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1225           outputContainer->Add(fhLambda0DispEta[i]);
1226           
1227           fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
1228                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1229           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1230           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1231           outputContainer->Add(fhLambda0DispPhi[i]);
1232           
1233         }
1234       }
1235     }
1236
1237     fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1238                            nptbins,ptmin,ptmax,20,0,20);
1239     fhNLocMaxPt ->SetYTitle("N maxima");
1240     fhNLocMaxPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1241     outputContainer->Add(fhNLocMaxPt) ;
1242
1243     for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1244     {
1245       fhNLocMaxPtSM[iSM] = new TH2F(Form("hNLocMaxPt_SM%d",iSM),Form("Number of local maxima in cluster, selected clusters in SM %d",iSM),
1246                                nptbins,ptmin,ptmax,20,0,20);
1247       fhNLocMaxPtSM[iSM] ->SetYTitle("N maxima");
1248       fhNLocMaxPtSM[iSM] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1249       outputContainer->Add(fhNLocMaxPtSM[iSM]) ;
1250     }
1251     
1252     if(fAnaType == kSSCalo)
1253     {
1254
1255       fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
1256                              nptbins,ptmin,ptmax,20,0,20);
1257       fhNLocMaxPtReject ->SetYTitle("N maxima");
1258       fhNLocMaxPtReject ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1259       outputContainer->Add(fhNLocMaxPtReject) ;
1260     }
1261     
1262     for (Int_t i = 0; i < 3; i++)
1263     {
1264       fhPtLambda0LocMax[i]  = new TH2F(Form("hPtLambda0LocMax%d",i+1),
1265                                       Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
1266                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1267       fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1268       fhPtLambda0LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1269       outputContainer->Add(fhPtLambda0LocMax[i]) ;
1270
1271       if(IsDataMC())
1272       {
1273         for(Int_t ipart = 0; ipart < 6; ipart++)
1274         {
1275           fhMCPtLambda0LocMax[ipart][i]  = new TH2F
1276           (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
1277            Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),pname[ipart].Data()),
1278            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1279           fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
1280           fhMCPtLambda0LocMax[ipart][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1281           outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
1282         }
1283       }
1284       
1285       fhPtLambda1LocMax[i]  = new TH2F(Form("hPtLambda1LocMax%d",i+1),
1286                                       Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}, %s",nlm[i].Data()),
1287                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1288       fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1289       fhPtLambda1LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1290       outputContainer->Add(fhPtLambda1LocMax[i]) ;
1291       
1292       fhPtDispersionLocMax[i]  = new TH2F(Form("hPtDispersionLocMax%d",i+1),
1293                                          Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion^{2}, %s",nlm[i].Data()),
1294                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1295       fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1296       fhPtDispersionLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1297       outputContainer->Add(fhPtDispersionLocMax[i]) ;
1298       
1299       if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1300       {
1301         fhPtDispEtaLocMax[i]  = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
1302                                         Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1303                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1304         fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1305         fhPtDispEtaLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1306         outputContainer->Add(fhPtDispEtaLocMax[i]) ;
1307         
1308         fhPtDispPhiLocMax[i]  = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
1309                                         Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1310                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1311         fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1312         fhPtDispPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1313         outputContainer->Add(fhPtDispPhiLocMax[i]) ;
1314         
1315         fhPtSumEtaPhiLocMax[i]  = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
1316                                           Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1317                                           nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1318         fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1319         fhPtSumEtaPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1320         outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
1321         
1322         fhPtDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
1323                                                Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1324                                                nptbins,ptmin,ptmax,200, -10,10);
1325         fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1326         fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1327         outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
1328         
1329         fhPtSphericityLocMax[i]  = new TH2F(Form("hPtSphericityLocMax%d",i+1),
1330                                            Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1331                                            nptbins,ptmin,ptmax,200, -1,1);
1332         fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1333         fhPtSphericityLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1334         outputContainer->Add(fhPtSphericityLocMax[i]) ;
1335       }
1336       
1337     }
1338     
1339     fhPtNCells  = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1340     fhPtNCells->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1341     fhPtNCells->SetYTitle("# of cells in cluster");
1342     outputContainer->Add(fhPtNCells);
1343     
1344     fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1345     fhPtTime->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1346     fhPtTime->SetYTitle("t (ns)");
1347     outputContainer->Add(fhPtTime);
1348     
1349   }
1350   
1351   
1352   fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1353   fhEPairDiffTime->SetXTitle("#it{E}_{pair} (GeV)");
1354   fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1355   outputContainer->Add(fhEPairDiffTime);
1356   
1357   if(fAnaType == kIMCalo)
1358   {
1359     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1360     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1361       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1362       "2 Local Maxima paired with more than 2 Local Maxima",
1363       "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1364     
1365     for (Int_t i = 0; i < 8 ; i++)
1366     {
1367       
1368       if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1369       
1370       fhMassPairLocMax[i]  = new TH2F
1371       (Form("MassPairLocMax%s",combiName[i].Data()),
1372        Form("#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1373        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1374       fhMassPairLocMax[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1375       fhMassPairLocMax[i]->SetXTitle("#it{E}_{pair} (GeV)");
1376       outputContainer->Add(fhMassPairLocMax[i]) ;
1377     }
1378   }
1379   
1380   if(fFillTMHisto)
1381   {
1382     fhTrackMatchedDEta  = new TH2F
1383     ("hTrackMatchedDEta",
1384      "d#eta of cluster-track vs cluster #it{p}_{T}",
1385      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1386     fhTrackMatchedDEta->SetYTitle("d#eta");
1387     fhTrackMatchedDEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1388     
1389     fhTrackMatchedDPhi  = new TH2F
1390     ("hTrackMatchedDPhi",
1391      "d#phi of cluster-track vs cluster #it{p}_{T}",
1392      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1393     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1394     fhTrackMatchedDPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1395     
1396     fhTrackMatchedDEtaDPhi  = new TH2F
1397     ("hTrackMatchedDEtaDPhi",
1398      "d#eta vs d#phi of cluster-track",
1399      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1400     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1401     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1402     
1403     outputContainer->Add(fhTrackMatchedDEta) ;
1404     outputContainer->Add(fhTrackMatchedDPhi) ;
1405     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1406
1407     fhTrackMatchedDEtaPos  = new TH2F
1408     ("hTrackMatchedDEtaPos",
1409      "d#eta of cluster-track vs cluster #it{p}_{T}",
1410      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1411     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1412     fhTrackMatchedDEtaPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1413     
1414     fhTrackMatchedDPhiPos  = new TH2F
1415     ("hTrackMatchedDPhiPos",
1416      "d#phi of cluster-track vs cluster #it{p}_{T}",
1417      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1418     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1419     fhTrackMatchedDPhiPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1420     
1421     fhTrackMatchedDEtaDPhiPos  = new TH2F
1422     ("hTrackMatchedDEtaDPhiPos",
1423      "d#eta vs d#phi of cluster-track",
1424      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1425     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1426     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1427     
1428     outputContainer->Add(fhTrackMatchedDEtaPos) ;
1429     outputContainer->Add(fhTrackMatchedDPhiPos) ;
1430     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1431
1432     fhTrackMatchedDEtaNeg  = new TH2F
1433     ("hTrackMatchedDEtaNeg",
1434      "d#eta of cluster-track vs cluster #it{p}_{T}",
1435      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1436     fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1437     fhTrackMatchedDEtaNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1438     
1439     fhTrackMatchedDPhiNeg  = new TH2F
1440     ("hTrackMatchedDPhiNeg",
1441      "d#phi of cluster-track vs cluster #it{p}_{T}",
1442      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1443     fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1444     fhTrackMatchedDPhiNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1445     
1446     fhTrackMatchedDEtaDPhiNeg  = new TH2F
1447     ("hTrackMatchedDEtaDPhiNeg",
1448      "d#eta vs d#phi of cluster-track",
1449      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1450     fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1451     fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1452     
1453     outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1454     outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1455     outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1456     
1457     fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1458     fhdEdx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1459     fhdEdx->SetYTitle("<#it{dE}/#it{dx}>");
1460     outputContainer->Add(fhdEdx);
1461     
1462     fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1463     fhEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1464     fhEOverP->SetYTitle("#it{E}/#it{p}");
1465     outputContainer->Add(fhEOverP);
1466     
1467     if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >=0)
1468     {
1469       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1470       fhEOverPNoTRD->SetXTitle("#it{E} (GeV)");
1471       fhEOverPNoTRD->SetYTitle("#it{E}/#it{p}");
1472       outputContainer->Add(fhEOverPNoTRD);
1473     }
1474     
1475     if(IsDataMC() && fFillTMHisto)
1476     {
1477       fhTrackMatchedMCParticlePt  = new TH2F
1478       ("hTrackMatchedMCParticlePt",
1479        "Origin of particle vs energy",
1480        nptbins,ptmin,ptmax,8,0,8);
1481       fhTrackMatchedMCParticlePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1482       //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
1483       
1484       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
1485       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
1486       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1487       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
1488       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1489       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1490       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1491       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1492       
1493       outputContainer->Add(fhTrackMatchedMCParticlePt);
1494       
1495       fhTrackMatchedMCParticleDEta  = new TH2F
1496       ("hTrackMatchedMCParticleDEta",
1497        "Origin of particle vs #eta residual",
1498        nresetabins,resetamin,resetamax,8,0,8);
1499       fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1500       //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1501       
1502       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1503       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1504       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1505       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1506       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1507       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1508       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1509       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1510       
1511       outputContainer->Add(fhTrackMatchedMCParticleDEta);
1512       
1513       fhTrackMatchedMCParticleDPhi  = new TH2F
1514       ("hTrackMatchedMCParticleDPhi",
1515        "Origin of particle vs #phi residual",
1516        nresphibins,resphimin,resphimax,8,0,8);
1517       fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1518       //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1519       
1520       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1521       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1522       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1523       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1524       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1525       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1526       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1527       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1528       
1529       outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1530       
1531       
1532     }
1533   }
1534   
1535   if(fFillWeightHistograms)
1536   {
1537     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1538                                      nptbins,ptmin,ptmax, 100,0,1.);
1539     fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1540     fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
1541     outputContainer->Add(fhECellClusterRatio);
1542     
1543     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1544                                         nptbins,ptmin,ptmax, 100,-10,0);
1545     fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1546     fhECellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1547     outputContainer->Add(fhECellClusterLogRatio);
1548     
1549     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1550                                         nptbins,ptmin,ptmax, 100,0,1.);
1551     fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1552     fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
1553     outputContainer->Add(fhEMaxCellClusterRatio);
1554     
1555     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1556                                            nptbins,ptmin,ptmax, 100,-10,0);
1557     fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1558     fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1559     outputContainer->Add(fhEMaxCellClusterLogRatio);
1560     
1561     for(Int_t iw = 0; iw < 14; iw++)
1562     {
1563       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
1564                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1565       fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1566       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1567       outputContainer->Add(fhLambda0ForW0[iw]);
1568       
1569       //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
1570       //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1571       //      fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1572       //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1573       //      outputContainer->Add(fhLambda1ForW0[iw]);
1574       
1575     }
1576   }
1577   
1578   if(IsDataMC())
1579   {
1580     // Origin
1581     
1582     fhMCPi0PtOrigin     = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1583     fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1584     fhMCPi0PtOrigin->SetYTitle("Origin");
1585     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1586     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1587     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1588     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1589     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1590     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1591     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1592     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1593     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1594     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1595     outputContainer->Add(fhMCPi0PtOrigin) ;
1596     
1597     fhMCEtaPtOrigin     = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1598     fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1599     fhMCEtaPtOrigin->SetYTitle("Origin");
1600     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1601     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1602     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1603     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1604     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1605     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1606     outputContainer->Add(fhMCEtaPtOrigin) ;
1607     
1608     fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1609     fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1610     fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1611     outputContainer->Add(fhMCPi0ProdVertex) ;
1612     
1613     fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1614     fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1615     fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1616     outputContainer->Add(fhMCEtaProdVertex) ;
1617     
1618     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1619     {
1620       fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #pi^{0} / E reco",
1621                                           nptbins,ptmin,ptmax,200,0,2);
1622       fhMCPi0PtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1623       fhMCPi0PtGenRecoFraction->SetYTitle("#it{E}^{#pi^{0} mother} / #it{E}^{rec}");
1624       outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1625       
1626       fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),#it{p}_{T} versus E primary #eta / E reco",
1627                                           nptbins,ptmin,ptmax,200,0,2);
1628       fhMCEtaPtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1629       fhMCEtaPtGenRecoFraction->SetYTitle("#it{E}^{ #eta mother} / #it{E}^{rec}");
1630       outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1631       
1632       fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1633       fhMCPi0DecayPt->SetYTitle("#it{N}");
1634       fhMCPi0DecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1635       outputContainer->Add(fhMCPi0DecayPt) ;
1636       
1637       fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta), #it{p}_{T} versus E primary  #gamma / #it{E} primary #pi^{0}",
1638                                         nptbins,ptmin,ptmax,100,0,1);
1639       fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/#it{c})");
1640       fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1641       outputContainer->Add(fhMCPi0DecayPtFraction) ;
1642       
1643       fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1644       fhMCEtaDecayPt->SetYTitle("#it{N}");
1645       fhMCEtaDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1646       outputContainer->Add(fhMCEtaDecayPt) ;
1647       
1648       fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), #it{p}_{T} versus E primary  #gamma / E primary #eta",
1649                                         nptbins,ptmin,ptmax,100,0,1);
1650       fhMCEtaDecayPtFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1651       fhMCEtaDecayPtFraction->SetYTitle("#it{E}^{gen} / #it{E}^{gen-mother}");
1652       outputContainer->Add(fhMCEtaDecayPtFraction) ;
1653       
1654       fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0})  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1655       fhMCOtherDecayPt->SetYTitle("#it{N}");
1656       fhMCOtherDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1657       outputContainer->Add(fhMCOtherDecayPt) ;
1658       
1659     }
1660     
1661     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1662        GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1663     {
1664       
1665       fhAnglePairMCPi0  = new TH2F
1666       ("AnglePairMCPi0",
1667        "Angle between decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1668       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1669       fhAnglePairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1670       outputContainer->Add(fhAnglePairMCPi0) ;
1671       
1672       if (fAnaType!= kSSCalo)
1673       {
1674         fhAnglePairMCEta  = new TH2F
1675         ("AnglePairMCEta",
1676          "Angle between decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1677         fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1678         fhAnglePairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1679         outputContainer->Add(fhAnglePairMCEta) ;
1680         
1681         fhMassPairMCPi0  = new TH2F
1682         ("MassPairMCPi0",
1683          "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1684         fhMassPairMCPi0->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1685         fhMassPairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1686         outputContainer->Add(fhMassPairMCPi0) ;
1687         
1688         fhMassPairMCEta  = new TH2F
1689         ("MassPairMCEta",
1690          "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1691         fhMassPairMCEta->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1692         fhMassPairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1693         outputContainer->Add(fhMassPairMCEta) ;
1694       }
1695       
1696       for(Int_t i = 0; i < 6; i++)
1697       {
1698         
1699         fhMCE[i]  = new TH1F
1700         (Form("hE_MC%s",pname[i].Data()),
1701          Form("Identified as #pi^{0} (#eta), cluster from %s",
1702               ptype[i].Data()),
1703          nptbins,ptmin,ptmax);
1704         fhMCE[i]->SetYTitle("#it{N}");
1705         fhMCE[i]->SetXTitle("#it{E} (GeV)");
1706         outputContainer->Add(fhMCE[i]) ;
1707         
1708         fhMCPt[i]  = new TH1F
1709         (Form("hPt_MC%s",pname[i].Data()),
1710          Form("Identified as #pi^{0} (#eta), cluster from %s",
1711               ptype[i].Data()),
1712          nptbins,ptmin,ptmax);
1713         fhMCPt[i]->SetYTitle("#it{N}");
1714         fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1715         outputContainer->Add(fhMCPt[i]) ;
1716         
1717         fhMCPtCentrality[i]  = new TH2F
1718         (Form("hPtCentrality_MC%s",pname[i].Data()),
1719          Form("Identified as #pi^{0} (#eta), cluster from %s",
1720               ptype[i].Data()),
1721          nptbins,ptmin,ptmax, 100,0,100);
1722         fhMCPtCentrality[i]->SetYTitle("centrality");
1723         fhMCPtCentrality[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1724         outputContainer->Add(fhMCPtCentrality[i]) ;
1725         
1726         if(fAnaType == kSSCalo)
1727         {
1728           fhMCNLocMaxPt[i] = new TH2F
1729           (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1730            Form("cluster from %s, #it{p}_{T} of cluster vs NLM, accepted",ptype[i].Data()),
1731            nptbins,ptmin,ptmax,20,0,20);
1732           fhMCNLocMaxPt[i] ->SetYTitle("#it{NLM}");
1733           fhMCNLocMaxPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1734           outputContainer->Add(fhMCNLocMaxPt[i]) ;
1735  
1736           fhMCNLocMaxPtReject[i] = new TH2F
1737           (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1738            Form("cluster from %s, #it{p}_{T} of cluster vs NLM, rejected",ptype[i].Data()),
1739            nptbins,ptmin,ptmax,20,0,20);
1740           fhMCNLocMaxPtReject[i] ->SetYTitle("#it{NLM}");
1741           fhMCNLocMaxPtReject[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1742           outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1743           
1744           fhMCEReject[i]  = new TH1F
1745           (Form("hEReject_MC%s",pname[i].Data()),
1746            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1747                 ptype[i].Data()),
1748            nptbins,ptmin,ptmax);
1749           fhMCEReject[i]->SetYTitle("#it{N}");
1750           fhMCEReject[i]->SetXTitle("#it{E} (GeV)");
1751           outputContainer->Add(fhMCEReject[i]) ;
1752           
1753           fhMCPtReject[i]  = new TH1F
1754           (Form("hPtReject_MC%s",pname[i].Data()),
1755            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1756                 ptype[i].Data()),
1757            nptbins,ptmin,ptmax);
1758           fhMCPtReject[i]->SetYTitle("#it{N}");
1759           fhMCPtReject[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1760           outputContainer->Add(fhMCPtReject[i]) ;
1761         }
1762         
1763         fhMCPtPhi[i]  = new TH2F
1764         (Form("hPtPhi_MC%s",pname[i].Data()),
1765          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1766          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1767         fhMCPtPhi[i]->SetYTitle("#phi");
1768         fhMCPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1769         outputContainer->Add(fhMCPtPhi[i]) ;
1770         
1771         fhMCPtEta[i]  = new TH2F
1772         (Form("hPtEta_MC%s",pname[i].Data()),
1773          Form("Identified as #pi^{0} (#eta), cluster from %s",
1774               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1775         fhMCPtEta[i]->SetYTitle("#eta");
1776         fhMCPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1777         outputContainer->Add(fhMCPtEta[i]) ;
1778         
1779         fhMCMassPt[i]  = new TH2F
1780         (Form("hMassPt_MC%s",pname[i].Data()),
1781          Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1782          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1783         fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1784         fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1785         outputContainer->Add(fhMCMassPt[i]) ;
1786         
1787         fhMCSelectedMassPt[i]  = new TH2F
1788         (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1789          Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1790          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1791         fhMCSelectedMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1792         fhMCSelectedMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1793         outputContainer->Add(fhMCSelectedMassPt[i]) ;
1794         
1795         if(fAnaType == kSSCalo)
1796         {
1797           fhMCMassPtNoOverlap[i]  = new TH2F
1798           (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1799            Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1800            nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1801           fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1802           fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1803           outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1804           
1805           fhMCSelectedMassPtNoOverlap[i]  = new TH2F
1806           (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1807            Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1808            nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1809           fhMCSelectedMassPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1810           fhMCSelectedMassPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1811           outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1812         }
1813         
1814         if( fFillSelectClHisto )
1815         {
1816           fhMCPtLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1817                                       Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
1818                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1819           fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1820           fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1821           outputContainer->Add(fhMCPtLambda0[i]) ;
1822           
1823           fhMCPtLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1824                                       Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
1825                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1826           fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1827           fhMCPtLambda1[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1828           outputContainer->Add(fhMCPtLambda1[i]) ;
1829           
1830           fhMCPtDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1831                                          Form("Selected pair, cluster from %s : #it{p}_{T} vs dispersion^{2}",ptype[i].Data()),
1832                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1833           fhMCPtDispersion[i]->SetYTitle("#it{D}^{2}");
1834           fhMCPtDispersion[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1835           outputContainer->Add(fhMCPtDispersion[i]) ;
1836           
1837           if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >= 0)
1838           {
1839             fhMCPtLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1840                                              Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1841                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1842             fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1843             fhMCPtLambda0NoTRD[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1844             outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
1845             
1846             if(!fFillOnlySimpleSSHisto)
1847             {
1848               fhMCPtDispEta[i]  = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
1849                                            Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
1850                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1851               fhMCPtDispEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1852               fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1853               outputContainer->Add(fhMCPtDispEta[i]);
1854               
1855               fhMCPtDispPhi[i]  = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
1856                                            Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
1857                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1858               fhMCPtDispPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1859               fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1860               outputContainer->Add(fhMCPtDispPhi[i]);
1861               
1862               fhMCPtSumEtaPhi[i]  = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
1863                                              Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",ptype[i].Data()),
1864                                              nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1865               fhMCPtSumEtaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1866               fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1867               outputContainer->Add(fhMCPtSumEtaPhi[i]);
1868               
1869               fhMCPtDispEtaPhiDiff[i]  = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
1870                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",ptype[i].Data()),
1871                                                   nptbins,ptmin,ptmax,200,-10,10);
1872               fhMCPtDispEtaPhiDiff[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1873               fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1874               outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
1875               
1876               fhMCPtSphericity[i]  = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
1877                                               Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
1878                                               nptbins,ptmin,ptmax, 200,-1,1);
1879               fhMCPtSphericity[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1880               fhMCPtSphericity[i]->SetYTitle("#it{s} = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1881               outputContainer->Add(fhMCPtSphericity[i]);
1882               
1883               for(Int_t ie = 0; ie < 7; ie++)
1884               {
1885                 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1886                                                       Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1887                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1888                 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1889                 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1890                 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1891                 
1892                 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1893                                                       Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1894                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1895                 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1896                 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1897                 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1898                 
1899                 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1900                                                       Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1901                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1902                 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1903                 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1904                 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1905                 
1906               }
1907             }
1908           }
1909           
1910           fhMCPtLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1911                                                     Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1912                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1913           fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1914           fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("#it{E} (GeV)");
1915           outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
1916           
1917           fhMCPtFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1918                                           Form("Selected pair, cluster from %s : #it{p}_{T} vs Max cell fraction of energy",ptype[i].Data()),
1919                                           nptbins,ptmin,ptmax,100,0,1);
1920           fhMCPtFracMaxCell[i]->SetYTitle("#it{Fraction}");
1921           fhMCPtFracMaxCell[i]->SetXTitle("#it{E} (GeV)");
1922           outputContainer->Add(fhMCPtFracMaxCell[i]) ;
1923           
1924         }//
1925       } // shower shape histo
1926       
1927     } //Not MC reader
1928   }//Histos with MC
1929   
1930   if(fAnaType==kSSCalo)
1931   {
1932     fhAsymmetry  = new TH2F ("hAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
1933                              nptbins,ptmin,ptmax, 200, -1,1);
1934     fhAsymmetry->SetXTitle("#it{E} (GeV)");
1935     fhAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
1936     outputContainer->Add(fhAsymmetry);
1937     
1938     fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
1939                                      nptbins,ptmin,ptmax, 200, -1,1);
1940     fhSelectedAsymmetry->SetXTitle("#it{E} (GeV)");
1941     fhSelectedAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
1942     outputContainer->Add(fhSelectedAsymmetry);
1943     
1944     fhSplitE  = new TH1F
1945     ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1946     fhSplitE->SetYTitle("counts");
1947     fhSplitE->SetXTitle("#it{E} (GeV)");
1948     outputContainer->Add(fhSplitE) ;
1949     
1950     fhSplitPt  = new TH1F
1951     ("hSplitPt","Selected #pi^{0} (#eta) pairs #it{p}_{T} sum of split sub-clusters",nptbins,ptmin,ptmax);
1952     fhSplitPt->SetYTitle("counts");
1953     fhSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1954     outputContainer->Add(fhSplitPt) ;
1955     
1956     
1957     fhSplitPtPhi  = new TH2F
1958     ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1959     fhSplitPtPhi->SetYTitle("#phi (rad)");
1960     fhSplitPtPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1961     outputContainer->Add(fhSplitPtPhi) ;
1962     
1963     fhSplitPtEta  = new TH2F
1964     ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1965     fhSplitPtEta->SetYTitle("#eta");
1966     fhSplitPtEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1967     outputContainer->Add(fhSplitPtEta) ;
1968     
1969     
1970     fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1971                                 nptbins,ptmin,ptmax,20,0,20);
1972     fhNLocMaxSplitPt ->SetYTitle("#it{NLM}");
1973     fhNLocMaxSplitPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1974     outputContainer->Add(fhNLocMaxSplitPt) ;
1975     
1976     
1977     fhMassSplitPt  = new TH2F
1978     ("hMassSplitPt","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
1979      nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1980     fhMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1981     fhMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1982     outputContainer->Add(fhMassSplitPt) ;
1983     
1984     fhSelectedMassSplitPt  = new TH2F
1985     ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
1986      nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1987     fhSelectedMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1988     fhSelectedMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1989     outputContainer->Add(fhSelectedMassSplitPt) ;
1990     
1991     if(IsDataMC())
1992     {
1993       fhMassSplitPtNoOverlap  = new TH2F
1994       ("hMassSplitPtNoOverlap","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
1995        nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1996       fhMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1997       fhMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1998       outputContainer->Add(fhMassSplitPtNoOverlap) ;
1999       
2000       fhSelectedMassSplitPtNoOverlap  = new TH2F
2001       ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
2002        nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2003       fhSelectedMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2004       fhSelectedMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2005       outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
2006
2007       
2008       fhMCPi0PtRecoPtPrim  = new TH2F
2009       ("hMCPi0PtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2010        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2011       fhMCPi0PtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2012       fhMCPi0PtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2013       outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
2014       
2015       fhMCPi0PtRecoPtPrimNoOverlap  = new TH2F
2016       ("hMCPi0PtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2017        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2018       fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2019       fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2020       outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
2021       
2022       fhMCPi0SelectedPtRecoPtPrim  = new TH2F
2023       ("hMCPi0SelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2024        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2025       fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2026       fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2027       outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
2028       
2029       fhMCPi0SelectedPtRecoPtPrimNoOverlap  = new TH2F
2030       ("hMCPi0SelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2031        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2032       fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2033       fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2034       outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
2035
2036       
2037       fhMCPi0SplitPtRecoPtPrim  = new TH2F
2038       ("hMCPi0SplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2039        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2040       fhMCPi0SplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2041       fhMCPi0SplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2042       outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
2043       
2044       fhMCPi0SplitPtRecoPtPrimNoOverlap  = new TH2F
2045       ("hMCPi0SplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2046        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2047       fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2048       fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2049       outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
2050       
2051       fhMCPi0SelectedSplitPtRecoPtPrim  = new TH2F
2052       ("hMCPi0SelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2053        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2054       fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2055       fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2056       outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
2057       
2058       fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap  = new TH2F
2059       ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2060        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2061       fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2062       fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2063       outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
2064
2065       fhMCEtaPtRecoPtPrim  = new TH2F
2066       ("hMCEtaPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2067        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2068       fhMCEtaPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2069       fhMCEtaPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2070       outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
2071       
2072       fhMCEtaPtRecoPtPrimNoOverlap  = new TH2F
2073       ("hMCEtaPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2074        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2075       fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2076       fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2077       outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
2078       
2079       fhMCEtaSelectedPtRecoPtPrim  = new TH2F
2080       ("hMCEtaSelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2081        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2082       fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2083       fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2084       outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
2085       
2086       fhMCEtaSelectedPtRecoPtPrimNoOverlap  = new TH2F
2087       ("hMCEtaSelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2088        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2089       fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2090       fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2091       outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
2092       
2093       
2094       fhMCEtaSplitPtRecoPtPrim  = new TH2F
2095       ("hMCEtaSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2096        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2097       fhMCEtaSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2098       fhMCEtaSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2099       outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
2100       
2101       fhMCEtaSplitPtRecoPtPrimNoOverlap  = new TH2F
2102       ("hMCEtaSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2103        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2104       fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2105       fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2106       outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
2107       
2108       fhMCEtaSelectedSplitPtRecoPtPrim  = new TH2F
2109       ("hMCEtaSelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2110        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2111       fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2112       fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2113       outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
2114       
2115       fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap  = new TH2F
2116       ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2117        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2118       fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2119       fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2120       outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
2121       
2122       
2123       for(Int_t inlm = 0; inlm < 3; inlm++)
2124       {
2125         fhMCPi0PtRecoPtPrimLocMax[inlm]  = new TH2F
2126         (Form("hMCPi0PtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2127          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2128         fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2129         fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2130         outputContainer->Add(fhMCPi0PtRecoPtPrimLocMax[inlm] ) ;
2131         
2132         fhMCPi0SelectedPtRecoPtPrimLocMax[inlm]  = new TH2F
2133         (Form("hMCPi0SelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2134          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2135         fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2136         fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2137         outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ) ;
2138         
2139         fhMCPi0SplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2140         (Form("hMCPi0SplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2141          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2142         fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2143         fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2144         outputContainer->Add(fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ) ;
2145         
2146         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2147         (Form("hMCPi0SelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2148          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2149         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2150         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2151         outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2152         
2153         fhMCEtaPtRecoPtPrimLocMax[inlm]  = new TH2F
2154         (Form("hMCEtaPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2155          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2156         fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2157         fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2158         outputContainer->Add(fhMCEtaPtRecoPtPrimLocMax[inlm] ) ;
2159         
2160         fhMCEtaSelectedPtRecoPtPrimLocMax[inlm]  = new TH2F
2161         (Form("hMCEtaSelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2162          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2163         fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2164         fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2165         outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ) ;
2166         
2167         fhMCEtaSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2168         (Form("hMCEtaSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2169          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2170         fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2171         fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2172         outputContainer->Add(fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ) ;
2173         
2174         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2175         (Form("hMCEtaSelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2176          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2177         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2178         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2179         outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2180         
2181       }
2182       
2183       for(Int_t i = 0; i< 6; i++)
2184       {
2185         fhMCPtAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
2186                                        Form("cluster from %s : #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",ptype[i].Data()),
2187                                        nptbins,ptmin,ptmax, 200,-1,1);
2188         fhMCPtAsymmetry[i]->SetXTitle("#it{E} (GeV)");
2189         fhMCPtAsymmetry[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2190         outputContainer->Add(fhMCPtAsymmetry[i]);
2191         
2192         fhMCSplitE[i]  = new TH1F
2193         (Form("hSplitE_MC%s",pname[i].Data()),
2194          Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2195          nptbins,ptmin,ptmax);
2196         fhMCSplitE[i]->SetYTitle("counts");
2197         fhMCSplitE[i]->SetXTitle("#it{E} (GeV)");
2198         outputContainer->Add(fhMCSplitE[i]) ;
2199         
2200         fhMCSplitPt[i]  = new TH1F
2201         (Form("hSplitPt_MC%s",pname[i].Data()),
2202          Form("cluster from %s, #it{p}_{T} sum of split sub-clusters",ptype[i].Data()),
2203          nptbins,ptmin,ptmax);
2204         fhMCSplitPt[i]->SetYTitle("counts");
2205         fhMCSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2206         outputContainer->Add(fhMCSplitPt[i]) ;
2207         
2208         
2209         fhMCSplitPtPhi[i]  = new TH2F
2210         (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2211          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2212          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2213         fhMCSplitPtPhi[i]->SetYTitle("#phi");
2214         fhMCSplitPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2215         outputContainer->Add(fhMCSplitPtPhi[i]) ;
2216         
2217         fhMCSplitPtEta[i]  = new TH2F
2218         (Form("hSplitPtEta_MC%s",pname[i].Data()),
2219          Form("Identified as #pi^{0} (#eta), cluster from %s",
2220               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2221         fhMCSplitPtEta[i]->SetYTitle("#eta");
2222         fhMCSplitPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2223         outputContainer->Add(fhMCSplitPtEta[i]) ;
2224         
2225         
2226         fhMCNLocMaxSplitPt[i] = new TH2F
2227         (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2228          Form("cluster from %s, #it{p}_{T} sum of split sub-clusters, for NLM",ptype[i].Data()),
2229          nptbins,ptmin,ptmax,20,0,20);
2230         fhMCNLocMaxSplitPt[i] ->SetYTitle("#it{NLM}");
2231         fhMCNLocMaxSplitPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2232         outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2233         
2234         fhMCMassSplitPt[i]  = new TH2F
2235         (Form("hMassSplitPt_MC%s",pname[i].Data()),
2236          Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2237          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2238         fhMCMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2239         fhMCMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2240         outputContainer->Add(fhMCMassSplitPt[i]) ;
2241         
2242         fhMCSelectedMassSplitPt[i]  = new TH2F
2243         (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2244          Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2245          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2246         fhMCSelectedMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2247         fhMCSelectedMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2248         outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2249
2250         fhMCMassSplitPtNoOverlap[i]  = new TH2F
2251         (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2252          Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2253          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2254         fhMCMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2255         fhMCMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2256         outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2257         
2258         fhMCSelectedMassSplitPtNoOverlap[i]  = new TH2F
2259         (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2260          Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2261          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2262         fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2263         fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2264         outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2265       }
2266     }
2267   }
2268   
2269   if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2270   {
2271     
2272     
2273     for(Int_t i = 0; i< 3; i++)
2274     {
2275       fhPtAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2276                                         Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ), %s",nlm[i].Data()),
2277                                         nptbins,ptmin,ptmax,200, -1,1);
2278       fhPtAsymmetryLocMax[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2279       fhPtAsymmetryLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2280       outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
2281     }
2282     
2283     for(Int_t ie = 0; ie< 7; ie++)
2284     {
2285       
2286       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2287                                          Form("#lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2288                                          ssbins,ssmin,ssmax , 200,-1,1);
2289       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2290       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2291       outputContainer->Add(fhAsymmetryLambda0[ie]);
2292       
2293       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2294                                          Form("#sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2295                                          ssbins,ssmin,ssmax , 200,-1,1);
2296       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2297       fhAsymmetryDispEta[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2298       outputContainer->Add(fhAsymmetryDispEta[ie]);
2299       
2300       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2301                                          Form("#sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2302                                          ssbins,ssmin,ssmax , 200,-1,1);
2303       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2304       fhAsymmetryDispPhi[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2305       outputContainer->Add(fhAsymmetryDispPhi[ie]);
2306     }
2307     
2308     
2309     if(IsDataMC())
2310     {
2311       for(Int_t i = 0; i< 6; i++)
2312       {
2313         for(Int_t ie = 0; ie < 7; ie++)
2314         {
2315           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2316                                                   Form("cluster from %s : #lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2317                                                   ssbins,ssmin,ssmax , 200,-1,1);
2318           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2319           fhMCAsymmetryLambda0[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2320           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2321           
2322           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2323                                                   Form("cluster from %s : #sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2324                                                   ssbins,ssmin,ssmax , 200,-1,1);
2325           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2326           fhMCAsymmetryDispEta[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2327           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2328           
2329           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2330                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2331                                                   ssbins,ssmin,ssmax , 200,-1,1);
2332           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2333           fhMCAsymmetryDispPhi[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2334           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2335         }
2336       }
2337     }
2338   }
2339   
2340   if(fFillPileUpHistograms)
2341   {
2342     
2343     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2344     
2345     for(Int_t i = 0 ; i < 7 ; i++)
2346     {
2347       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2348                                    Form("Selected #pi^{0} (#eta) #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2349       fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2350       outputContainer->Add(fhPtPileUp[i]);
2351       
2352       fhPtCellTimePileUp[i]  = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2353                                              Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2354                                              nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
2355       fhPtCellTimePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2356       fhPtCellTimePileUp[i]->SetYTitle("#it{t}_{cell} (ns)");
2357       outputContainer->Add(fhPtCellTimePileUp[i]);
2358       
2359       fhPtTimeDiffPileUp[i]  = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2360                                              Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2361                                              nptbins,ptmin,ptmax,400,-200,200);
2362       fhPtTimeDiffPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2363       fhPtTimeDiffPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
2364       outputContainer->Add(fhPtTimeDiffPileUp[i]);
2365
2366     }
2367     
2368     fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","#it{t} of cluster vs #it{E} of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2369     fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2370     fhTimePtNoCut->SetYTitle("#it{t} (ns)");
2371     outputContainer->Add(fhTimePtNoCut);
2372     
2373     fhTimePtSPD  = new TH2F ("hTimePt_SPD","#it{t} of cluster vs #it{E} of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2374     fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2375     fhTimePtSPD->SetYTitle("#it{t} (ns)");
2376     outputContainer->Add(fhTimePtSPD);
2377     
2378     fhTimePtSPDMulti  = new TH2F ("hTimePt_SPDMulti","time of cluster vs #it{E} of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2379     fhTimePtSPDMulti->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2380     fhTimePtSPDMulti->SetYTitle("#it{t} (ns)");
2381     outputContainer->Add(fhTimePtSPDMulti);
2382     
2383     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","#it{t} of cluster vs #it{N} pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2384     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2385     fhTimeNPileUpVertSPD->SetXTitle("#it{t} (ns)");
2386     outputContainer->Add(fhTimeNPileUpVertSPD);
2387     
2388     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","#it{t} of cluster vs #it{N} pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
2389     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2390     fhTimeNPileUpVertTrack->SetXTitle("#it{t} (ns)");
2391     outputContainer->Add(fhTimeNPileUpVertTrack);
2392     
2393     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","#it{t} of cluster vs #it{N} constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2394     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2395     fhTimeNPileUpVertContributors->SetXTitle("#it{t} (ns)");
2396     outputContainer->Add(fhTimeNPileUpVertContributors);
2397     
2398     fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - main SPD vertex",ntimptbins,timemin,timemax,100,0,50);
2399     fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{Z} (cm) ");
2400     fhTimePileUpMainVertexZDistance->SetXTitle("#it{t} (ns)");
2401     outputContainer->Add(fhTimePileUpMainVertexZDistance);
2402     
2403     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - z diamond",ntimptbins,timemin,timemax,100,0,50);
2404     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{Z} (cm) ");
2405     fhTimePileUpMainVertexZDiamond->SetXTitle("#it{t} (ns)");
2406     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2407                 
2408                 fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","#it{p}_{T} of cluster vs #it{N} pile-up SPD vertex",
2409                                                                                                                                          nptbins,ptmin,ptmax,20,0,20);
2410                 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2411                 fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2412                 outputContainer->Add(fhPtNPileUpSPDVtx);
2413           
2414                 fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","#it{p}_{T} of cluster vs #it{N} pile-up Tracks vertex",
2415                                                                                                                                          nptbins,ptmin,ptmax, 20,0,20 );
2416                 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2417                 fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2418                 outputContainer->Add(fhPtNPileUpTrkVtx);
2419                 
2420                 fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","#it{p}_{T} of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2421                                           nptbins,ptmin,ptmax,20,0,20);
2422                 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2423                 fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2424                 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2425           
2426                 fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","#it{p}_{T} of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2427                                                                                                                                                                         nptbins,ptmin,ptmax, 20,0,20 );
2428                 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2429                 fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2430                 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2431     
2432     fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","#it{p}_{T} of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2433                                            nptbins,ptmin,ptmax,20,0,20);
2434                 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2435                 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2436                 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2437           
2438                 fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","#it{p}_{T} of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2439                                            nptbins,ptmin,ptmax, 20,0,20 );
2440                 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2441                 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2442                 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2443     
2444   }
2445   
2446   //Keep neutral meson selection histograms if requiered
2447   //Setting done in AliNeutralMesonSelection
2448   
2449   if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2450   {
2451     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2452     
2453     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2454       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2455     
2456     delete nmsHistos;
2457   }
2458   
2459   return outputContainer ;
2460   
2461 }
2462
2463 //_____________________________________________
2464 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2465 {
2466   
2467   // Assign mc index depending on MC bit set
2468   
2469   if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
2470   {
2471     return kmcPi0 ;
2472   }//pi0
2473   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
2474   {
2475     return kmcEta ;
2476   }//eta
2477   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2478              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2479   {
2480     return kmcConversion ;
2481   }//conversion photon
2482   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
2483   {
2484     return kmcPhoton ;
2485   }//photon   no conversion
2486   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2487   {
2488     return kmcElectron ;
2489   }//electron
2490   else
2491   {
2492     return kmcHadron ;
2493   }//other particles
2494   
2495 }
2496
2497 //__________________________________________________________________
2498 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2499                                        AliAODPWG4Particle * photon2,
2500                                        Int_t & label, Int_t & tag)
2501 {
2502   // Check the labels of pare in case mother was same pi0 or eta
2503   // Set the new AOD accordingly
2504   
2505   Int_t  label1 = photon1->GetLabel();
2506   Int_t  label2 = photon2->GetLabel();
2507   
2508   if(label1 < 0 || label2 < 0 ) return ;
2509   
2510   //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2511   //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2512   Int_t tag1 = photon1->GetTag();
2513   Int_t tag2 = photon2->GetTag();
2514   
2515   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2516   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2517        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
2518       (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2519        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
2520      )
2521   {
2522     
2523     //Check if pi0/eta mother is the same
2524     if(GetReader()->ReadStack())
2525     {
2526       if(label1>=0)
2527       {
2528         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2529         label1 = mother1->GetFirstMother();
2530         //mother1 = GetMCStack()->Particle(label1);//pi0
2531       }
2532       if(label2>=0)
2533       {
2534         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2535         label2 = mother2->GetFirstMother();
2536         //mother2 = GetMCStack()->Particle(label2);//pi0
2537       }
2538     } // STACK
2539     else if(GetReader()->ReadAODMCParticles())
2540     {//&& (input > -1)){
2541       if(label1>=0)
2542       {
2543         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2544         label1 = mother1->GetMother();
2545         //mother1 = GetMCStack()->Particle(label1);//pi0
2546       }
2547       if(label2>=0)
2548       {
2549         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2550         label2 = mother2->GetMother();
2551         //mother2 = GetMCStack()->Particle(label2);//pi0
2552       }
2553     }// AOD
2554     
2555     //printf("mother1 %d, mother2 %d\n",label1,label2);
2556     if( label1 == label2 && label1>=0 )
2557     {
2558       
2559       label = label1;
2560       
2561       TLorentzVector mom1 = *(photon1->Momentum());
2562       TLorentzVector mom2 = *(photon2->Momentum());
2563       
2564       Double_t angle = mom2.Angle(mom1.Vect());
2565       Double_t mass  = (mom1+mom2).M();
2566       Double_t epair = (mom1+mom2).E();
2567       
2568       if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2569       {
2570         fhMassPairMCPi0 ->Fill(epair,mass);
2571         fhAnglePairMCPi0->Fill(epair,angle);
2572         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2573       }
2574       else
2575       {
2576         fhMassPairMCEta ->Fill(epair,mass);
2577         fhAnglePairMCEta->Fill(epair,angle);
2578         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2579       }
2580       
2581     } // same label
2582   } // both from eta or pi0 decay
2583   
2584 }
2585
2586 //____________________________________________________________________________
2587 void AliAnaPi0EbE::Init()
2588 {
2589   //Init
2590   //Do some checks
2591   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2592     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2593     abort();
2594   }
2595   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2596     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2597     abort();
2598   }
2599   
2600 }
2601
2602 //____________________________________________________________________________
2603 void AliAnaPi0EbE::InitParameters()
2604 {
2605   //Initialize the parameters of the analysis.
2606   AddToHistogramsName("AnaPi0EbE_");
2607   
2608   fInputAODGammaConvName = "PhotonsCTS" ;
2609   fAnaType = kIMCalo ;
2610   fCalorimeter = "EMCAL" ;
2611   fMinDist  = 2.;
2612   fMinDist2 = 4.;
2613   fMinDist3 = 5.;
2614   
2615   fNLMECutMin[0] = 10.;
2616   fNLMECutMin[1] = 6. ;
2617   fNLMECutMin[2] = 6. ;
2618 }
2619
2620 //__________________________________________________________________
2621 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
2622 {
2623   //Do analysis and fill aods
2624   
2625   switch(fAnaType)
2626   {
2627     case kIMCalo:
2628       MakeInvMassInCalorimeter();
2629       break;
2630       
2631     case kSSCalo:
2632       MakeShowerShapeIdentification();
2633       break;
2634       
2635     case kIMCaloTracks:
2636       MakeInvMassInCalorimeterAndCTS();
2637       break;
2638       
2639   }
2640 }
2641
2642 //____________________________________________
2643 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
2644 {
2645   //Do analysis and fill aods
2646   //Search for the photon decay in calorimeters
2647   //Read photon list from AOD, produced in class AliAnaPhoton
2648   //Check if 2 photons have the mass of the pi0.
2649   
2650   TLorentzVector mom1;
2651   TLorentzVector mom2;
2652   TLorentzVector mom ;
2653   
2654   Int_t tag   = 0;
2655   Int_t label = 0;
2656   
2657   if(!GetInputAODBranch()){
2658     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2659     abort();
2660   }
2661   
2662   //Get shower shape information of clusters
2663   TObjArray *clusters = 0;
2664   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2665   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2666   
2667   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2668     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2669     
2670     //Vertex cut in case of mixed events
2671     Int_t evtIndex1 = 0 ;
2672     if(GetMixedEvent())
2673       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2674     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
2675     mom1 = *(photon1->Momentum());
2676     
2677     //Get original cluster, to recover some information
2678     Int_t iclus = -1;
2679     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2680     
2681     if(!cluster1){
2682       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2683       return;
2684     }
2685     
2686     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2687     {
2688       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2689       
2690       Int_t evtIndex2 = 0 ;
2691       if(GetMixedEvent())
2692         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2693       
2694       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2695         continue ;
2696       
2697       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
2698       
2699       mom2 = *(photon2->Momentum());
2700       
2701       //Get original cluster, to recover some information
2702       Int_t iclus2 = -1;
2703       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2704       // start new loop from iclus1+1 to gain some time
2705
2706       if(!cluster2)
2707       {
2708         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2709         return;
2710       }
2711       
2712       Float_t e1    = photon1->E();
2713       Float_t e2    = photon2->E();
2714       
2715       //Select clusters with good time window difference
2716       Float_t tof1  = cluster1->GetTOF()*1e9;;
2717       Float_t tof2  = cluster2->GetTOF()*1e9;;
2718       Double_t t12diff = tof1-tof2;
2719       fhEPairDiffTime->Fill(e1+e2,    t12diff);
2720       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2721       
2722       //Play with the MC stack if available
2723       if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2724       
2725       // Check the invariant mass for different selection on the local maxima
2726       // Name of AOD method TO BE FIXED
2727       Int_t nMaxima1 = photon1->GetFiducialArea();
2728       Int_t nMaxima2 = photon2->GetFiducialArea();
2729       
2730       mom = mom1+mom2;
2731       
2732       Double_t mass  = mom.M();
2733       Double_t epair = mom.E();
2734       
2735       if(nMaxima1==nMaxima2)
2736       {
2737         if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2738         else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2739         else                 fhMassPairLocMax[2]->Fill(epair,mass);
2740       }
2741       else if(nMaxima1==1 || nMaxima2==1)
2742       {
2743         if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2744         else                             fhMassPairLocMax[4]->Fill(epair,mass);
2745       }
2746       else
2747         fhMassPairLocMax[5]->Fill(epair,mass);
2748       
2749       // combinations with SS axis cut and NLM cut
2750       if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2751       if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2752       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2753       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2754       
2755       //Skip events with too few or too many  NLM
2756       if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2757       
2758       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2759       
2760       //Mass of all pairs
2761       fhMass->Fill(epair,mass);
2762       
2763       //Select good pair (good phi, pt cuts, aperture and invariant mass)
2764       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2765       {
2766         if(GetDebug()>1)
2767           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",
2768                  mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
2769         
2770         //Fill some histograms about shower shape
2771         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2772         {
2773           FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2774           FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2775         }
2776         
2777         // Tag both photons as decay, !!careful!! since in case of SideBand analysis, also they will be tagged.
2778         photon1->SetTagged(kTRUE);
2779         photon2->SetTagged(kTRUE);
2780         
2781         fhPtDecay->Fill(photon1->Pt());
2782         fhEDecay ->Fill(photon1->E() );
2783         
2784         fhPtDecay->Fill(photon2->Pt());
2785         fhEDecay ->Fill(photon2->E() );
2786         
2787         //Mass of selected pairs
2788         fhSelectedMass->Fill(epair,mass);
2789         
2790         // Fill histograms to undertand pile-up before other cuts applied
2791         // Remember to relax time cuts in the reader
2792         FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2793         
2794         //Create AOD for analysis
2795         
2796         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2797         
2798         if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2799         else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2800         else
2801         {
2802           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2803           return ;
2804         }
2805         pi0.SetDetector(photon1->GetDetector());
2806         
2807         // MC
2808         pi0.SetLabel(label);
2809         pi0.SetTag(tag);
2810         
2811         //Set the indeces of the original caloclusters
2812         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2813         //pi0.SetInputFileIndex(input);
2814         
2815         AddAODParticle(pi0);
2816         
2817       }//pi0
2818       
2819     }//2n photon loop
2820     
2821   }//1st photon loop
2822   
2823   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2824   
2825 }
2826
2827 //__________________________________________________
2828 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2829 {
2830   //Do analysis and fill aods
2831   //Search for the photon decay in calorimeters
2832   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2833   //Check if 2 photons have the mass of the pi0.
2834   
2835   TLorentzVector mom1;
2836   TLorentzVector mom2;
2837   TLorentzVector mom ;
2838   Int_t tag   = 0;
2839   Int_t label = 0;
2840   Int_t evtIndex = 0;
2841   
2842   // Check calorimeter input
2843   if(!GetInputAODBranch())
2844   {
2845     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2846     abort();
2847   }
2848   
2849   // Get the array with conversion photons
2850   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2851   if(!inputAODGammaConv)
2852   {
2853     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2854     
2855     if(!inputAODGammaConv)
2856     {
2857       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2858       
2859       return;
2860     }
2861   }
2862   
2863   //Get shower shape information of clusters
2864   TObjArray *clusters = 0;
2865   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2866   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2867   
2868   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
2869   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2870   if(nCTS<=0 || nCalo <=0)
2871   {
2872     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2873     return;
2874   }
2875   
2876   if(GetDebug() > 1)
2877     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2878   
2879   // Do the loop, first calo, second CTS
2880   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++)
2881   {
2882     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2883     mom1 = *(photon1->Momentum());
2884     
2885     //Get original cluster, to recover some information
2886     Int_t iclus = -1;
2887     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2888     
2889     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++)
2890     {
2891       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2892       
2893       if(GetMixedEvent())
2894         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2895       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
2896       
2897       mom2 = *(photon2->Momentum());
2898       
2899       mom = mom1+mom2;
2900       
2901       Double_t mass  = mom.M();
2902       Double_t epair = mom.E();
2903       
2904       Int_t nMaxima = photon1->GetFiducialArea();
2905       if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2906       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2907       else                fhMassPairLocMax[2]->Fill(epair,mass);
2908       
2909       if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2910       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2911       
2912       //Play with the MC stack if available
2913       if(IsDataMC())
2914       {
2915         Int_t   label2 = photon2->GetLabel();
2916         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2917         
2918         HasPairSameMCMother(photon1, photon2, label, tag) ;
2919       }
2920       
2921       //Mass of selected pairs
2922       fhMass->Fill(epair,mass);
2923       
2924       //Select good pair (good phi, pt cuts, aperture and invariant mass)
2925       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2926       {
2927         if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",
2928                                   mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
2929         
2930         //Fill some histograms about shower shape
2931         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2932         {
2933           FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
2934         }
2935         
2936         // Tag both photons as decay, !!careful!! since in case of SideBand analysis, also they will be tagged.
2937         photon1->SetTagged(kTRUE);
2938         photon2->SetTagged(kTRUE);
2939         
2940         fhPtDecay->Fill(photon1->Pt());
2941         fhEDecay ->Fill(photon1->E() );
2942         
2943         //Mass of selected pairs
2944         fhSelectedMass->Fill(epair,mass);
2945         
2946         // Fill histograms to undertand pile-up before other cuts applied
2947         // Remember to relax time cuts in the reader
2948         if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2949
2950         //Create AOD for analysis
2951         
2952         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2953         
2954         if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2955         else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2956         else
2957         {
2958           printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2959           return ;
2960         }
2961         pi0.SetDetector(photon1->GetDetector());
2962         
2963         // MC
2964         pi0.SetLabel(label);
2965         pi0.SetTag(tag);
2966         
2967         //Set the indeces of the original tracks or caloclusters
2968         pi0.SetCaloLabel (photon1->GetCaloLabel(0) , -1);
2969         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2970         //pi0.SetInputFileIndex(input);
2971         
2972         AddAODParticle(pi0);
2973         
2974       }//pi0
2975     }//2n photon loop
2976     
2977   }//1st photon loop
2978   
2979   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2980   
2981 }
2982
2983
2984 //_________________________________________________
2985 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
2986 {
2987   //Search for pi0 in fCalorimeter with shower shape analysis
2988   
2989   TObjArray * pl        = 0x0;
2990   AliVCaloCells * cells = 0x0;
2991   //Select the Calorimeter of the photon
2992   if      (fCalorimeter == "PHOS" )
2993   {
2994     pl    = GetPHOSClusters();
2995     cells = GetPHOSCells();
2996   }
2997   else if (fCalorimeter == "EMCAL")
2998   {
2999     pl    = GetEMCALClusters();
3000     cells = GetEMCALCells();
3001   }
3002   
3003   if(!pl)
3004   {
3005     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
3006     return;
3007   }
3008         
3009   TLorentzVector mom ;
3010   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
3011   {
3012     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
3013     
3014     Int_t evtIndex = 0 ;
3015     if (GetMixedEvent())
3016     {
3017       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3018     }
3019     
3020     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
3021     
3022     //Get Momentum vector,
3023     Double_t vertex[]={0,0,0};
3024     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3025     {
3026       calo->GetMomentum(mom,GetVertex(evtIndex)) ;