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