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