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