coverity, check if pointer was null
[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     fUseSplitAsyCut(kFALSE),
55     fTimeCutMin(-10000),           fTimeCutMax(10000),
56     fFillPileUpHistograms(0),
57     fFillWeightHistograms(kFALSE), fFillTMHisto(0),              
58     fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),
59     fInputAODGammaConvName(""),
60     // Histograms
61     fhPt(0),                       fhE(0),                    
62     fhEEta(0),                     fhEPhi(0),                    fhEtaPhi(0),
63     fhPtReject(0),                 fhEReject(0),                    
64     fhEEtaReject(0),               fhEPhiReject(0),              fhEtaPhiReject(0),
65     fhMass(0),                     fhAsymmetry(0), 
66     fhSelectedMass(0),             fhSelectedAsymmetry(0),
67     fhPtDecay(0),                  fhEDecay(0),  
68     // Shower shape histos
69     fhEDispersion(0),              fhELambda0(0),                fhELambda1(0), 
70     fhELambda0NoTRD(0),            fhELambda0FracMaxCellCut(0),  
71     fhEFracMaxCell(0),             fhEFracMaxCellNoTRD(0),            
72     fhENCells(0),                  fhETime(0),                   fhEPairDiffTime(0),
73     fhDispEtaE(0),                 fhDispPhiE(0),
74     fhSumEtaE(0),                  fhSumPhiE(0),                 fhSumEtaPhiE(0),
75     fhDispEtaPhiDiffE(0),          fhSphericityE(0),           
76
77     // MC histos
78     fhMCE(),                       fhMCPt(),        
79     fhMCPhi(),                     fhMCEta(),
80     fhMCEReject(),                 fhMCPtReject(),       
81     fhMCPi0PtGenRecoFraction(0),   fhMCEtaPtGenRecoFraction(0),
82     fhMCPi0DecayPt(0),             fhMCPi0DecayPtFraction(0),      
83     fhMCEtaDecayPt(0),             fhMCEtaDecayPtFraction(0),
84     fhMCOtherDecayPt(0),           
85     fhMassPairMCPi0(0),            fhMassPairMCEta(0),
86     fhAnglePairMCPi0(0),           fhAnglePairMCEta(0),
87     // Weight studies
88     fhECellClusterRatio(0),        fhECellClusterLogRatio(0),                 
89     fhEMaxCellClusterRatio(0),     fhEMaxCellClusterLogRatio(0),
90     fhTrackMatchedDEta(0),         fhTrackMatchedDPhi(0),        fhTrackMatchedDEtaDPhi(0),
91     fhTrackMatchedMCParticle(0),   fhdEdx(0),                     
92     fhEOverP(0),                   fhEOverPNoTRD(0),                
93     // Number of local maxima in cluster
94     fhNLocMax(0),
95     // PileUp
96     fhTimeENoCut(0),                    fhTimeESPD(0),           fhTimeESPDMulti(0),
97     fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
98     fhTimeNPileUpVertContributors(0),
99     fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
100 {
101   //default ctor
102   
103   for(Int_t i = 0; i < 6; i++)
104   {
105     fhMCE              [i] = 0;
106     fhMCPt             [i] = 0;
107     fhMCPhi            [i] = 0;                   
108     fhMCEta            [i] = 0;
109     fhEMCLambda0       [i] = 0;
110     fhEMCLambda0NoTRD  [i] = 0;
111     fhEMCLambda0FracMaxCellCut[i]= 0;
112     fhEMCFracMaxCell   [i] = 0;
113     fhEMCLambda1       [i] = 0;
114     fhEMCDispersion    [i] = 0;
115     
116     fhMCEDispEta       [i] = 0;
117     fhMCEDispPhi       [i] = 0;
118     fhMCESumEtaPhi     [i] = 0;
119     fhMCEDispEtaPhiDiff[i] = 0;
120     fhMCESphericity    [i] = 0;    
121     fhMCEAsymmetry     [i] = 0;          
122
123     for(Int_t j = 0; j < 7; j++)
124     {    
125       fhMCLambda0DispEta    [j][i] = 0;
126       fhMCLambda0DispPhi    [j][i] = 0;
127       fhMCDispEtaDispPhi    [j][i] = 0; 
128       fhMCAsymmetryLambda0  [j][i] = 0;    
129       fhMCAsymmetryDispEta  [j][i] = 0; 
130       fhMCAsymmetryDispPhi  [j][i] = 0;
131     }
132   }
133   
134   for(Int_t j = 0; j < 7; j++)
135   {  
136     fhLambda0DispEta    [j] = 0;
137     fhLambda0DispPhi    [j] = 0;
138     fhDispEtaDispPhi    [j] = 0; 
139     fhAsymmetryLambda0  [j] = 0;    
140     fhAsymmetryDispEta  [j] = 0; 
141     fhAsymmetryDispPhi  [j] = 0;
142     
143     fhPtPi0PileUp       [j] = 0;
144   }
145   
146   for(Int_t i = 0; i < 3; i++)
147   {
148     fhELambda0LocMax       [i] = 0;
149     fhELambda1LocMax       [i] = 0;
150     fhEDispersionLocMax    [i] = 0;  
151     fhEDispEtaLocMax       [i] = 0;  
152     fhEDispPhiLocMax       [i] = 0;  
153     fhESumEtaPhiLocMax     [i] = 0;
154     fhEDispEtaPhiDiffLocMax[i] = 0;
155     fhESphericityLocMax    [i] = 0;
156     fhEAsymmetryLocMax     [i] = 0;
157   }
158   
159   //Weight studies
160   for(Int_t i =0; i < 14; i++){
161     fhLambda0ForW0[i] = 0;
162     //fhLambda1ForW0[i] = 0;
163     if(i<8)fhMassPairLocMax[i] = 0;
164   }
165   
166   //Initialize parameters
167   InitParameters();
168   
169 }
170
171 //_______________________________________________________________________________
172 void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time) 
173 {
174   // Fill some histograms to understand pile-up
175   if(!fFillPileUpHistograms) return;
176   
177   //printf("E %f, time %f\n",energy,time);
178   AliVEvent * event = GetReader()->GetInputEvent();
179   
180   fhTimeENoCut->Fill(energy,time);
181   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
182   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
183   
184   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
185   
186   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
187   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
188   
189   // N pile up vertices
190   Int_t nVerticesSPD    = -1;
191   Int_t nVerticesTracks = -1;
192   
193   if      (esdEv)
194   {
195     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
196     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
197     
198   }//ESD
199   else if (aodEv)
200   {
201     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
202     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
203   }//AOD
204   
205   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
206   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
207   
208   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n", 
209   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
210   
211   Int_t ncont = -1;
212   Float_t z1 = -1, z2 = -1;
213   Float_t diamZ = -1;
214   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
215   {
216     if      (esdEv)
217     {
218       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
219       ncont=pv->GetNContributors();
220       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
221       z2 = pv->GetZ();
222       diamZ = esdEv->GetDiamondZ();
223     }//ESD
224     else if (aodEv)
225     {
226       AliAODVertex *pv=aodEv->GetVertex(iVert);
227       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
228       ncont=pv->GetNContributors();
229       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
230       z2=pv->GetZ();
231       diamZ = aodEv->GetDiamondZ();
232     }// AOD
233     
234     Double_t distZ  = TMath::Abs(z2-z1);
235     diamZ  = TMath::Abs(z2-diamZ);
236     
237     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
238     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
239     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
240     
241   }// loop
242 }
243
244
245 //___________________________________________________________________________________________
246 void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
247 {
248  // Fill histograms that do not pass the identification (SS case only)
249   
250   Float_t ener  = mom.E();
251   Float_t pt    = mom.Pt();
252   Float_t phi   = mom.Phi();
253   if(phi < 0) phi+=TMath::TwoPi();
254   Float_t eta = mom.Eta();
255   
256   fhPtReject     ->Fill(pt);
257   fhEReject      ->Fill(ener);
258   
259   fhEEtaReject   ->Fill(ener,eta);
260   fhEPhiReject   ->Fill(ener,phi);
261   fhEtaPhiReject ->Fill(eta,phi);
262   
263   if(IsDataMC())
264   {
265     Int_t mcIndex = GetMCIndex(mctag);
266     fhMCEReject  [mcIndex] ->Fill(ener);
267     fhMCPtReject [mcIndex] ->Fill(pt);
268   }    
269 }
270
271 //_____________________________________________________________________________________
272 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, 
273                                                  const Int_t nMaxima,
274                                                  const Int_t tag, 
275                                                  const Float_t asy)
276 {
277   // Fill shower shape, timing and other histograms for selected clusters from decay
278   
279   Float_t e    = cluster->E();
280   Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
281   Float_t l0   = cluster->GetM02();
282   Float_t l1   = cluster->GetM20(); 
283   Int_t   nSM  = GetModuleNumber(cluster);
284
285   Int_t ebin = -1;
286   if      (e < 2 ) ebin = 0;
287   else if (e < 4 ) ebin = 1;
288   else if (e < 6 ) ebin = 2;
289   else if (e < 10) ebin = 3;
290   else if (e < 15) ebin = 4;  
291   else if (e < 20) ebin = 5;  
292   else             ebin = 6;  
293
294   Int_t indexMax = -1;
295   if     (nMaxima==1) indexMax = 0 ;
296   else if(nMaxima==2) indexMax = 1 ; 
297   else                indexMax = 2 ; 
298   
299   
300   AliVCaloCells * cell = 0x0; 
301   if(fCalorimeter == "PHOS") 
302     cell = GetPHOSCells();
303   else                        
304     cell = GetEMCALCells();
305   
306   Float_t maxCellFraction = 0;
307   GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
308   fhEFracMaxCell->Fill(e,maxCellFraction);  
309   
310   FillWeightHistograms(cluster);
311   
312   fhEDispersion->Fill(e, disp);   
313   fhELambda0   ->Fill(e, l0  );  
314   fhELambda1   ->Fill(e, l1  );  
315   
316   Float_t ll0  = 0., ll1  = 0.;
317   Float_t dispp= 0., dEta = 0., dPhi    = 0.; 
318   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
319   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
320   {
321     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
322                                                                                  ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
323     
324     fhDispEtaE        -> Fill(e,dEta);
325     fhDispPhiE        -> Fill(e,dPhi);
326     fhSumEtaE         -> Fill(e,sEta);
327     fhSumPhiE         -> Fill(e,sPhi);
328     fhSumEtaPhiE      -> Fill(e,sEtaPhi);
329     fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
330     if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
331     
332     fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
333     fhLambda0DispEta[ebin]->Fill(l0  ,dEta);
334     fhLambda0DispPhi[ebin]->Fill(l0  ,dPhi);
335     
336     if (fAnaType==kSSCalo)
337     {
338       // Asymmetry histograms
339       fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
340       fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
341       fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
342     }
343   }  
344   
345   fhNLocMax->Fill(e,nMaxima);
346
347   fhELambda0LocMax   [indexMax]->Fill(e,l0); 
348   fhELambda1LocMax   [indexMax]->Fill(e,l1);
349   fhEDispersionLocMax[indexMax]->Fill(e,disp);
350   
351   if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto) 
352   {
353     fhEDispEtaLocMax       [indexMax]-> Fill(e,dEta);
354     fhEDispPhiLocMax       [indexMax]-> Fill(e,dPhi);
355     fhESumEtaPhiLocMax     [indexMax]-> Fill(e,sEtaPhi);
356     fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
357     if(dEta+dPhi>0)       fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
358     if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e  ,asy);
359     
360   }
361   
362   if(fCalorimeter=="EMCAL" && nSM < 6) 
363   {
364     fhELambda0NoTRD->Fill(e, l0  );
365     fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);  
366   }
367   
368   if(maxCellFraction < 0.5) 
369     fhELambda0FracMaxCellCut->Fill(e, l0  );  
370   
371   fhETime  ->Fill(e, cluster->GetTOF()*1.e9);
372   fhENCells->Fill(e, cluster->GetNCells());
373   
374   // Fill Track matching control histograms
375   if(fFillTMHisto)
376   {
377     Float_t dZ  = cluster->GetTrackDz();
378     Float_t dR  = cluster->GetTrackDx();
379
380     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
381     {
382       dR = 2000., dZ = 2000.;
383       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
384     }    
385     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
386
387     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
388     {
389       fhTrackMatchedDEta->Fill(e,dZ);
390       fhTrackMatchedDPhi->Fill(e,dR);
391       if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);      
392     }
393     
394     // Check dEdx and E/p of matched clusters
395     
396     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
397     {
398       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
399       
400       if(track) 
401       {
402         Float_t dEdx = track->GetTPCsignal();
403         fhdEdx->Fill(e, dEdx);
404         
405         Float_t eOverp = e/track->P();
406         fhEOverP->Fill(e,  eOverp);
407         
408         if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e,  eOverp);
409
410       }
411       //else 
412       //  printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
413       
414
415       
416       if(IsDataMC())
417       {
418         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
419         {
420           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
421                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(e, 2.5 );
422           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(e, 0.5 );
423           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(e, 1.5 );
424           else                                                                                 fhTrackMatchedMCParticle->Fill(e, 3.5 );
425           
426         }
427         else
428         {
429           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
430                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(e, 6.5 );
431           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(e, 4.5 );
432           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(e, 5.5 );
433           else                                                                                 fhTrackMatchedMCParticle->Fill(e, 7.5 );
434         }        
435       }  // MC              
436     }
437   }// Track matching histograms   
438   
439   if(IsDataMC()) 
440   {
441     Int_t mcIndex = GetMCIndex(tag);
442     
443     fhEMCLambda0[mcIndex]    ->Fill(e, l0);
444     fhEMCLambda1[mcIndex]    ->Fill(e, l1);
445     fhEMCDispersion[mcIndex] ->Fill(e, disp);
446     fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction); 
447     
448     if(fCalorimeter=="EMCAL" && nSM < 6) 
449       fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0  );
450     
451     if(maxCellFraction < 0.5) 
452       fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0  );  
453     
454     if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
455     {
456       fhMCEDispEta        [mcIndex]-> Fill(e,dEta);
457       fhMCEDispPhi        [mcIndex]-> Fill(e,dPhi);
458       fhMCESumEtaPhi      [mcIndex]-> Fill(e,sEtaPhi);
459       fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
460       if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));  
461
462       if (fAnaType==kSSCalo)
463       {
464         fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
465         fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
466         fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
467       }
468       
469       fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
470       fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0  ,dEta);
471       fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0  ,dPhi);
472       
473     }
474     
475   }//MC
476   
477 }
478
479 //________________________________________________________
480 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
481 {
482   // Calculate weights and fill histograms
483   
484   if(!fFillWeightHistograms || GetMixedEvent()) return;
485   
486   AliVCaloCells* cells = 0;
487   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
488   else                        cells = GetPHOSCells();
489   
490   // First recalculate energy in case non linearity was applied
491   Float_t  energy = 0;
492   Float_t  ampMax = 0;  
493   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) 
494   {
495     
496     Int_t id       = clus->GetCellsAbsId()[ipos];
497     
498     //Recalibrate cell energy if needed
499     Float_t amp = cells->GetCellAmplitude(id);
500     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
501     
502     energy    += amp;
503     
504     if(amp> ampMax) 
505       ampMax = amp;
506     
507   } // energy loop       
508   
509   if(energy <=0 ) 
510   {
511     printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
512     return;
513   }
514   
515   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
516   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
517   
518   //Get the ratio and log ratio to all cells in cluster
519   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) 
520   {
521     Int_t id       = clus->GetCellsAbsId()[ipos];
522     
523     //Recalibrate cell energy if needed
524     Float_t amp = cells->GetCellAmplitude(id);
525     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
526     
527     fhECellClusterRatio   ->Fill(energy,amp/energy);
528     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
529   }        
530   
531   //Recalculate shower shape for different W0
532   if(fCalorimeter=="EMCAL"){
533     
534     Float_t l0org = clus->GetM02();
535     Float_t l1org = clus->GetM20();
536     Float_t dorg  = clus->GetDispersion();
537     
538     for(Int_t iw = 0; iw < 14; iw++)
539     {
540       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5); 
541       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
542       
543       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
544       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
545       
546     } // w0 loop
547     
548     // Set the original values back
549     clus->SetM02(l0org);
550     clus->SetM20(l1org);
551     clus->SetDispersion(dorg);
552     
553   }// EMCAL
554 }
555
556 //__________________________________________
557 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
558 {       
559         //Save parameters used for analysis
560   TString parList ; //this will be list of parameters used for this analysis.
561   const Int_t buffersize = 255;
562   char onePar[buffersize] ;
563   
564   snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
565   parList+=onePar ;     
566   snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
567   parList+=onePar ;
568   
569   if(fAnaType == kSSCalo)
570   {
571     snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
572     parList+=onePar ;
573     snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
574     parList+=onePar ;
575     snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
576     parList+=onePar ;
577     snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
578     parList+=onePar ;
579   }
580   
581   //Get parameters set in base class.
582   parList += GetBaseParametersList() ;
583   
584   //Get parameters set in PID class.
585   if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
586   
587   return new TObjString(parList) ;
588 }
589
590 //_____________________________________________
591 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
592 {  
593   // Create histograms to be saved in output file and 
594   // store them in outputContainer
595   TList * outputContainer = new TList() ; 
596   outputContainer->SetName("Pi0EbEHistos") ; 
597   
598   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
599   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
600   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
601   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
602   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
603   Int_t tbins    = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin   = GetHistogramRanges()->GetHistoTimeMin();
604   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax   = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin   = GetHistogramRanges()->GetHistoNClusterCellMin(); 
605
606   Int_t   nmassbins   = GetHistogramRanges()->GetHistoMassBins();            
607   Float_t massmin     = GetHistogramRanges()->GetHistoMassMin();              
608   Float_t massmax     = GetHistogramRanges()->GetHistoMassMax();
609   
610   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();          
611   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();          
612   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
613   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();          
614   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();          
615   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
616   
617   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         
618   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();         
619   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
620   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       
621   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();       
622   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
623   
624   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();         
625   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();         
626   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();      
627   
628   TString nlm[]   ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
629   TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"}; 
630   TString pname[] ={"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};  
631   Int_t   bin[]   = {0,2,4,6,10,15,20,100}; // energy bins
632
633   fhPt  = new TH1F("hPt","Number of identified  #pi^{0} (#eta) decay",nptbins,ptmin,ptmax); 
634   fhPt->SetYTitle("N");
635   fhPt->SetXTitle("p_{T} (GeV/c)");
636   outputContainer->Add(fhPt) ; 
637   
638   fhE  = new TH1F("hE","Number of identified  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax); 
639   fhE->SetYTitle("N");
640   fhE->SetXTitle("E (GeV)");
641   outputContainer->Add(fhE) ; 
642   
643   fhEPhi  = new TH2F
644   ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
645   fhEPhi->SetYTitle("#phi (rad)");
646   fhEPhi->SetXTitle("E (GeV)");
647   outputContainer->Add(fhEPhi) ; 
648   
649   fhEEta  = new TH2F
650   ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
651   fhEEta->SetYTitle("#eta");
652   fhEEta->SetXTitle("E (GeV)");
653   outputContainer->Add(fhEEta) ; 
654   
655   fhEtaPhi  = new TH2F
656   ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax); 
657   fhEtaPhi->SetYTitle("#phi (rad)");
658   fhEtaPhi->SetXTitle("#eta");
659   outputContainer->Add(fhEtaPhi) ; 
660   
661   if(fAnaType == kSSCalo)
662   {
663     fhPtReject  = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax); 
664     fhPtReject->SetYTitle("N");
665     fhPtReject->SetXTitle("p_{T} (GeV/c)");
666     outputContainer->Add(fhPtReject) ; 
667     
668     fhEReject  = new TH1F("hEReject","Number of rejected as  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax); 
669     fhEReject->SetYTitle("N");
670     fhEReject->SetXTitle("E (GeV)");
671     outputContainer->Add(fhEReject) ; 
672     
673     fhEPhiReject  = new TH2F
674     ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
675     fhEPhiReject->SetYTitle("#phi (rad)");
676     fhEPhiReject->SetXTitle("E (GeV)");
677     outputContainer->Add(fhEPhiReject) ; 
678     
679     fhEEtaReject  = new TH2F
680     ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
681     fhEEtaReject->SetYTitle("#eta");
682     fhEEtaReject->SetXTitle("E (GeV)");
683     outputContainer->Add(fhEEtaReject) ; 
684     
685     fhEtaPhiReject  = new TH2F
686     ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax); 
687     fhEtaPhiReject->SetYTitle("#phi (rad)");
688     fhEtaPhiReject->SetXTitle("#eta");
689     outputContainer->Add(fhEtaPhiReject) ; 
690   }
691   
692   fhMass  = new TH2F
693   ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax); 
694   fhMass->SetYTitle("mass (GeV/c^{2})");
695   fhMass->SetXTitle("E (GeV)");
696   outputContainer->Add(fhMass) ; 
697   
698   fhSelectedMass  = new TH2F
699   ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax); 
700   fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
701   fhSelectedMass->SetXTitle("E (GeV)");
702   outputContainer->Add(fhSelectedMass) ; 
703   
704   if(fAnaType != kSSCalo)
705   {
706     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax); 
707     fhPtDecay->SetYTitle("N");
708     fhPtDecay->SetXTitle("p_{T} (GeV/c)");
709     outputContainer->Add(fhPtDecay) ; 
710     
711     fhEDecay  = new TH1F("hEDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax); 
712     fhEDecay->SetYTitle("N");
713     fhEDecay->SetXTitle("E (GeV)");
714     outputContainer->Add(fhEDecay) ;   
715   }
716   
717   ////////
718   
719   if( fFillSelectClHisto )
720   {
721     
722     fhEDispersion  = new TH2F
723     ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
724     fhEDispersion->SetYTitle("D^{2}");
725     fhEDispersion->SetXTitle("E (GeV)");
726     outputContainer->Add(fhEDispersion) ; 
727     
728     fhELambda0  = new TH2F
729     ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
730     fhELambda0->SetYTitle("#lambda_{0}^{2}");
731     fhELambda0->SetXTitle("E (GeV)");
732     outputContainer->Add(fhELambda0) ; 
733
734     fhELambda1  = new TH2F
735     ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
736     fhELambda1->SetYTitle("#lambda_{1}^{2}");
737     fhELambda1->SetXTitle("E (GeV)");
738     outputContainer->Add(fhELambda1) ; 
739         
740     fhELambda0FracMaxCellCut  = new TH2F
741     ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
742     fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
743     fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
744     outputContainer->Add(fhELambda0FracMaxCellCut) ; 
745
746     fhEFracMaxCell  = new TH2F
747     ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1); 
748     fhEFracMaxCell->SetYTitle("Fraction");
749     fhEFracMaxCell->SetXTitle("E (GeV)");
750     outputContainer->Add(fhEFracMaxCell) ; 
751     
752     if(fCalorimeter=="EMCAL")
753     {
754       fhELambda0NoTRD  = new TH2F
755       ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
756       fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
757       fhELambda0NoTRD->SetXTitle("E (GeV)");
758       outputContainer->Add(fhELambda0NoTRD) ; 
759       
760       fhEFracMaxCellNoTRD  = new TH2F
761       ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1); 
762       fhEFracMaxCellNoTRD->SetYTitle("Fraction");
763       fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
764       outputContainer->Add(fhEFracMaxCellNoTRD) ; 
765       
766       if(!fFillOnlySimpleSSHisto)
767       {
768         fhDispEtaE  = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
769         fhDispEtaE->SetXTitle("E (GeV)");
770         fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
771         outputContainer->Add(fhDispEtaE);     
772         
773         fhDispPhiE  = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
774         fhDispPhiE->SetXTitle("E (GeV)");
775         fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
776         outputContainer->Add(fhDispPhiE);  
777         
778         fhSumEtaE  = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
779         fhSumEtaE->SetXTitle("E (GeV)");
780         fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
781         outputContainer->Add(fhSumEtaE);     
782         
783         fhSumPhiE  = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",  
784                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
785         fhSumPhiE->SetXTitle("E (GeV)");
786         fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
787         outputContainer->Add(fhSumPhiE);  
788         
789         fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",  
790                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
791         fhSumEtaPhiE->SetXTitle("E (GeV)");
792         fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
793         outputContainer->Add(fhSumEtaPhiE);
794         
795         fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E", 
796                                        nptbins,ptmin,ptmax,200, -10,10); 
797         fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
798         fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
799         outputContainer->Add(fhDispEtaPhiDiffE);    
800         
801         fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",  
802                                    nptbins,ptmin,ptmax, 200, -1,1); 
803         fhSphericityE->SetXTitle("E (GeV)");
804         fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
805         outputContainer->Add(fhSphericityE);
806         
807         for(Int_t i = 0; i < 7; i++)
808         {
809           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]), 
810                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
811           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
812           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
813           outputContainer->Add(fhDispEtaDispPhi[i]); 
814           
815           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]), 
816                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
817           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
818           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
819           outputContainer->Add(fhLambda0DispEta[i]);       
820           
821           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]), 
822                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
823           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
824           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
825           outputContainer->Add(fhLambda0DispPhi[i]); 
826           
827         }
828       }
829     }    
830     
831     fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
832                          nptbins,ptmin,ptmax,10,0,10); 
833     fhNLocMax ->SetYTitle("N maxima");
834     fhNLocMax ->SetXTitle("E (GeV)");
835     outputContainer->Add(fhNLocMax) ;  
836     
837     for (Int_t i = 0; i < 3; i++) 
838     {
839       fhELambda0LocMax[i]  = new TH2F(Form("hELambda0LocMax%d",i+1),
840                                       Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
841                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
842       fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
843       fhELambda0LocMax[i]->SetXTitle("E (GeV)");
844       outputContainer->Add(fhELambda0LocMax[i]) ; 
845       
846       fhELambda1LocMax[i]  = new TH2F(Form("hELambda1LocMax%d",i+1),
847                                       Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
848                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
849       fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
850       fhELambda1LocMax[i]->SetXTitle("E (GeV)");
851       outputContainer->Add(fhELambda1LocMax[i]) ; 
852       
853       fhEDispersionLocMax[i]  = new TH2F(Form("hEDispersionLocMax%d",i+1),
854                                          Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
855                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
856       fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
857       fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
858       outputContainer->Add(fhEDispersionLocMax[i]) ; 
859       
860       if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
861       {
862         fhEDispEtaLocMax[i]  = new TH2F(Form("hEDispEtaLocMax%d",i+1),
863                                         Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
864                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
865         fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
866         fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
867         outputContainer->Add(fhEDispEtaLocMax[i]) ; 
868         
869         fhEDispPhiLocMax[i]  = new TH2F(Form("hEDispPhiLocMax%d",i+1),
870                                         Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
871                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
872         fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
873         fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
874         outputContainer->Add(fhEDispPhiLocMax[i]) ; 
875         
876         fhESumEtaPhiLocMax[i]  = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
877                                           Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
878                                           nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax); 
879         fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
880         fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
881         outputContainer->Add(fhESumEtaPhiLocMax[i]) ; 
882         
883         fhEDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
884                                                Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
885                                                nptbins,ptmin,ptmax,200, -10,10); 
886         fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
887         fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
888         outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ; 
889         
890         fhESphericityLocMax[i]  = new TH2F(Form("hESphericityLocMax%d",i+1),
891                                            Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
892                                            nptbins,ptmin,ptmax,200, -1,1); 
893         fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
894         fhESphericityLocMax[i]->SetXTitle("E (GeV)");
895         outputContainer->Add(fhESphericityLocMax[i]) ;
896       }
897        
898     }
899       
900     fhENCells  = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax); 
901     fhENCells->SetXTitle("E (GeV)");
902     fhENCells->SetYTitle("# of cells in cluster");
903     outputContainer->Add(fhENCells);  
904     
905     fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
906     fhETime->SetXTitle("E (GeV)");
907     fhETime->SetYTitle("t (ns)");
908     outputContainer->Add(fhETime);    
909     
910   }
911   
912   if(fAnaType == kIMCalo)
913   {
914     fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
915     fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
916     fhEPairDiffTime->SetYTitle("#Delta t (ns)");
917     outputContainer->Add(fhEPairDiffTime);
918     
919     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
920     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
921       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
922       "2 Local Maxima paired with more than 2 Local Maxima",
923       "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
924
925     for (Int_t i = 0; i < 8 ; i++) 
926     {
927
928       if (fAnaType == kIMCaloTracks && i > 2 ) continue ; 
929
930       fhMassPairLocMax[i]  = new TH2F
931       (Form("MassPairLocMax%s",combiName[i].Data()),
932        Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
933        nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
934       fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
935       fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
936       outputContainer->Add(fhMassPairLocMax[i]) ; 
937     }
938   }
939   
940   if(fFillTMHisto)
941   {
942     fhTrackMatchedDEta  = new TH2F
943     ("hTrackMatchedDEta",
944      "d#eta of cluster-track vs cluster energy",
945      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
946     fhTrackMatchedDEta->SetYTitle("d#eta");
947     fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
948     
949     fhTrackMatchedDPhi  = new TH2F
950     ("hTrackMatchedDPhi",
951      "d#phi of cluster-track vs cluster energy",
952      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
953     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
954     fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
955     
956     fhTrackMatchedDEtaDPhi  = new TH2F
957     ("hTrackMatchedDEtaDPhi",
958      "d#eta vs d#phi of cluster-track vs cluster energy",
959      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax); 
960     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
961     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");   
962     
963     outputContainer->Add(fhTrackMatchedDEta) ; 
964     outputContainer->Add(fhTrackMatchedDPhi) ;
965     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
966     
967     fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax); 
968     fhdEdx->SetXTitle("E (GeV)");
969     fhdEdx->SetYTitle("<dE/dx>");
970     outputContainer->Add(fhdEdx);  
971     
972     fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax); 
973     fhEOverP->SetXTitle("E (GeV)");
974     fhEOverP->SetYTitle("E/p");
975     outputContainer->Add(fhEOverP); 
976     
977     if(fCalorimeter=="EMCAL")
978     {
979       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax); 
980       fhEOverPNoTRD->SetXTitle("E (GeV)");
981       fhEOverPNoTRD->SetYTitle("E/p");
982       outputContainer->Add(fhEOverPNoTRD);   
983     }   
984     
985     if(IsDataMC() && fFillTMHisto)
986     {
987       fhTrackMatchedMCParticle  = new TH2F
988       ("hTrackMatchedMCParticle",
989        "Origin of particle vs energy",
990        nptbins,ptmin,ptmax,8,0,8); 
991       fhTrackMatchedMCParticle->SetXTitle("E (GeV)");   
992       //fhTrackMatchedMCParticle->SetYTitle("Particle type");
993       
994       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
995       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
996       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
997       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
998       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
999       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1000       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1001       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1002       
1003       outputContainer->Add(fhTrackMatchedMCParticle);   
1004     }
1005   }  
1006   
1007   if(fFillWeightHistograms)
1008   {
1009     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1010                                      nptbins,ptmin,ptmax, 100,0,1.); 
1011     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1012     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1013     outputContainer->Add(fhECellClusterRatio);
1014     
1015     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1016                                         nptbins,ptmin,ptmax, 100,-10,0); 
1017     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1018     fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1019     outputContainer->Add(fhECellClusterLogRatio);
1020     
1021     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1022                                         nptbins,ptmin,ptmax, 100,0,1.); 
1023     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1024     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1025     outputContainer->Add(fhEMaxCellClusterRatio);
1026     
1027     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1028                                            nptbins,ptmin,ptmax, 100,-10,0); 
1029     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1030     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1031     outputContainer->Add(fhEMaxCellClusterLogRatio);
1032     
1033     for(Int_t iw = 0; iw < 14; iw++)
1034     {
1035       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),
1036                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1037       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1038       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1039       outputContainer->Add(fhLambda0ForW0[iw]); 
1040       
1041 //      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),
1042 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1043 //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1044 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1045 //      outputContainer->Add(fhLambda1ForW0[iw]); 
1046       
1047     }
1048   }  
1049   
1050   if(IsDataMC()) 
1051   {
1052     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1053     {
1054       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",
1055                                    nptbins,ptmin,ptmax,200,0,2); 
1056       fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1057       fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1058       outputContainer->Add(fhMCPi0PtGenRecoFraction) ; 
1059             
1060       fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1061                                    nptbins,ptmin,ptmax,200,0,2); 
1062       fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1063       fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1064       outputContainer->Add(fhMCEtaPtGenRecoFraction) ; 
1065       
1066       fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax); 
1067       fhMCPi0DecayPt->SetYTitle("N");
1068       fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1069       outputContainer->Add(fhMCPi0DecayPt) ; 
1070       
1071       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}",
1072                                         nptbins,ptmin,ptmax,100,0,1); 
1073       fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1074       fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1075       outputContainer->Add(fhMCPi0DecayPtFraction) ; 
1076       
1077       fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax); 
1078       fhMCEtaDecayPt->SetYTitle("N");
1079       fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1080       outputContainer->Add(fhMCEtaDecayPt) ; 
1081       
1082       fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #eta",
1083                                         nptbins,ptmin,ptmax,100,0,1); 
1084       fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1085       fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1086       outputContainer->Add(fhMCEtaDecayPtFraction) ; 
1087       
1088       fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0})  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax); 
1089       fhMCOtherDecayPt->SetYTitle("N");
1090       fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1091       outputContainer->Add(fhMCOtherDecayPt) ; 
1092       
1093     }
1094        
1095     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
1096        GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1097     {
1098       
1099       fhAnglePairMCPi0  = new TH2F
1100       ("AnglePairMCPi0",
1101        "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5); 
1102       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1103       fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1104       outputContainer->Add(fhAnglePairMCPi0) ; 
1105
1106       if (fAnaType!= kSSCalo)
1107       {
1108         fhAnglePairMCEta  = new TH2F
1109         ("AnglePairMCEta",
1110          "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5); 
1111         fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1112         fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1113         outputContainer->Add(fhAnglePairMCEta) ; 
1114         
1115         fhMassPairMCPi0  = new TH2F
1116         ("MassPairMCPi0",
1117          "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
1118         fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1119         fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1120         outputContainer->Add(fhMassPairMCPi0) ; 
1121         
1122         fhMassPairMCEta  = new TH2F
1123         ("MassPairMCEta",
1124          "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
1125         fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1126         fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1127         outputContainer->Add(fhMassPairMCEta) ; 
1128       }
1129       
1130       for(Int_t i = 0; i < 6; i++)
1131       { 
1132         
1133         fhMCE[i]  = new TH1F
1134         (Form("hE_MC%s",pname[i].Data()),
1135          Form("Identified as #pi^{0} (#eta), cluster from %s",
1136               ptype[i].Data()),
1137          nptbins,ptmin,ptmax); 
1138         fhMCE[i]->SetYTitle("N");
1139         fhMCE[i]->SetXTitle("E (GeV)");
1140         outputContainer->Add(fhMCE[i]) ; 
1141         
1142         fhMCPt[i]  = new TH1F
1143         (Form("hPt_MC%s",pname[i].Data()),
1144          Form("Identified as #pi^{0} (#eta), cluster from %s",
1145               ptype[i].Data()),
1146          nptbins,ptmin,ptmax); 
1147         fhMCPt[i]->SetYTitle("N");
1148         fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1149         outputContainer->Add(fhMCPt[i]) ; 
1150         
1151         if(fAnaType == kSSCalo)
1152         {
1153           fhMCEReject[i]  = new TH1F
1154           (Form("hEReject_MC%s",pname[i].Data()),
1155            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1156                 ptype[i].Data()),
1157            nptbins,ptmin,ptmax); 
1158           fhMCEReject[i]->SetYTitle("N");
1159           fhMCEReject[i]->SetXTitle("E (GeV)");
1160           outputContainer->Add(fhMCEReject[i]) ; 
1161           
1162           fhMCPtReject[i]  = new TH1F
1163           (Form("hPtReject_MC%s",pname[i].Data()),
1164            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1165                 ptype[i].Data()),
1166            nptbins,ptmin,ptmax); 
1167           fhMCPtReject[i]->SetYTitle("N");
1168           fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1169           outputContainer->Add(fhMCPtReject[i]) ;           
1170         }
1171         
1172         fhMCPhi[i]  = new TH2F
1173         (Form("hPhi_MC%s",pname[i].Data()),
1174          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1175          nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
1176         fhMCPhi[i]->SetYTitle("#phi");
1177         fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1178         outputContainer->Add(fhMCPhi[i]) ; 
1179         
1180         fhMCEta[i]  = new TH2F
1181         (Form("hEta_MC%s",pname[i].Data()),
1182          Form("Identified as #pi^{0} (#eta), cluster from %s",
1183               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax); 
1184         fhMCEta[i]->SetYTitle("#eta");
1185         fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1186         outputContainer->Add(fhMCEta[i]) ;
1187         
1188         
1189         if( fFillSelectClHisto )
1190         {
1191           fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1192                                       Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1193                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1194           fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1195           fhEMCLambda0[i]->SetXTitle("E (GeV)");
1196           outputContainer->Add(fhEMCLambda0[i]) ; 
1197           
1198           fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1199                                       Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1200                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1201           fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1202           fhEMCLambda1[i]->SetXTitle("E (GeV)");
1203           outputContainer->Add(fhEMCLambda1[i]) ; 
1204           
1205           fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1206                                          Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1207                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1208           fhEMCDispersion[i]->SetYTitle("D^{2}");
1209           fhEMCDispersion[i]->SetXTitle("E (GeV)");
1210           outputContainer->Add(fhEMCDispersion[i]) ; 
1211           
1212           if(fCalorimeter=="EMCAL")
1213           {
1214             fhEMCLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1215                                              Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1216                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1217             fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1218             fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1219             outputContainer->Add(fhEMCLambda0NoTRD[i]) ; 
1220             
1221             if(!fFillOnlySimpleSSHisto)
1222             {
1223               fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1224                                            Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1225                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
1226               fhMCEDispEta[i]->SetXTitle("E (GeV)");
1227               fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1228               outputContainer->Add(fhMCEDispEta[i]);     
1229               
1230               fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1231                                            Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1232                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
1233               fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1234               fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1235               outputContainer->Add(fhMCEDispPhi[i]);  
1236               
1237               fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1238                                              Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),  
1239                                              nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
1240               fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1241               fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1242               outputContainer->Add(fhMCESumEtaPhi[i]);
1243               
1244               fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1245                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),  
1246                                                   nptbins,ptmin,ptmax,200,-10,10); 
1247               fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1248               fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1249               outputContainer->Add(fhMCEDispEtaPhiDiff[i]);    
1250               
1251               fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1252                                               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()),  
1253                                               nptbins,ptmin,ptmax, 200,-1,1); 
1254               fhMCESphericity[i]->SetXTitle("E (GeV)");
1255               fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1256               outputContainer->Add(fhMCESphericity[i]);
1257               
1258               for(Int_t ie = 0; ie < 7; ie++)
1259               {
1260                 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1261                                                       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]), 
1262                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
1263                 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1264                 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1265                 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]); 
1266                 
1267                 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1268                                                       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]), 
1269                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
1270                 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1271                 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1272                 outputContainer->Add(fhMCLambda0DispEta[ie][i]);       
1273                 
1274                 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1275                                                       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]), 
1276                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
1277                 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1278                 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1279                 outputContainer->Add(fhMCLambda0DispPhi[ie][i]); 
1280                 
1281               }            
1282             }
1283           }
1284           
1285           fhEMCLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1286                                                     Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1287                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1288           fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1289           fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1290           outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ; 
1291           
1292           fhEMCFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1293                                           Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1294                                           nptbins,ptmin,ptmax,100,0,1); 
1295           fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1296           fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1297           outputContainer->Add(fhEMCFracMaxCell[i]) ;           
1298           
1299         }//
1300       } // shower shape histo
1301       
1302     } //Not MC reader
1303   }//Histos with MC
1304   
1305   if(fAnaType==kSSCalo)
1306   {
1307     fhAsymmetry  = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",  
1308                               nptbins,ptmin,ptmax, 200, -1,1); 
1309     fhAsymmetry->SetXTitle("E (GeV)");
1310     fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1311     outputContainer->Add(fhAsymmetry);
1312     
1313     fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",  
1314                              nptbins,ptmin,ptmax, 200, -1,1); 
1315     fhSelectedAsymmetry->SetXTitle("E (GeV)");
1316     fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1317     outputContainer->Add(fhSelectedAsymmetry);
1318     
1319     if(IsDataMC())
1320     {
1321       for(Int_t i = 0; i< 6; i++)
1322       {
1323         fhMCEAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1324                                        Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),  
1325                                        nptbins,ptmin,ptmax, 200,-1,1); 
1326         fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1327         fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1328         outputContainer->Add(fhMCEAsymmetry[i]);
1329       } 
1330     }
1331   }
1332   
1333   if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
1334   {
1335     
1336     
1337     for(Int_t i = 0; i< 3; i++)
1338     {
1339       fhEAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1340                                         Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1341                                         nptbins,ptmin,ptmax,200, -1,1); 
1342       fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1343       fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1344       outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1345     }
1346     
1347     for(Int_t ie = 0; ie< 7; ie++)
1348     {
1349       
1350       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1351                                          Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]), 
1352                                          ssbins,ssmin,ssmax , 200,-1,1); 
1353       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1354       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1355       outputContainer->Add(fhAsymmetryLambda0[ie]); 
1356       
1357       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1358                                          Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]), 
1359                                          ssbins,ssmin,ssmax , 200,-1,1); 
1360       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1361       fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1362       outputContainer->Add(fhAsymmetryDispEta[ie]); 
1363       
1364       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1365                                          Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]), 
1366                                          ssbins,ssmin,ssmax , 200,-1,1); 
1367       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1368       fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1369       outputContainer->Add(fhAsymmetryDispPhi[ie]);           
1370     }        
1371     
1372     
1373     if(IsDataMC()) 
1374     {
1375       for(Int_t i = 0; i< 6; i++)
1376       {
1377         for(Int_t ie = 0; ie < 7; ie++)
1378         {
1379           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1380                                                   Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]), 
1381                                                   ssbins,ssmin,ssmax , 200,-1,1); 
1382           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1383           fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1384           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]); 
1385           
1386           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1387                                                   Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]), 
1388                                                   ssbins,ssmin,ssmax , 200,-1,1); 
1389           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1390           fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1391           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]); 
1392           
1393           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1394                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]), 
1395                                                   ssbins,ssmin,ssmax , 200,-1,1); 
1396           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1397           fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1398           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);     
1399         }        
1400       }
1401     }
1402   }
1403   
1404   if(fFillPileUpHistograms)
1405   {
1406     
1407     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1408
1409     for(Int_t i = 0 ; i < 7 ; i++)
1410     {
1411       fhPtPi0PileUp[i]  = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
1412                                       Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1413       fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
1414       outputContainer->Add(fhPtPi0PileUp[i]);
1415     }
1416     
1417     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1418     fhTimeENoCut->SetXTitle("E (GeV)");
1419     fhTimeENoCut->SetYTitle("time (ns)");
1420     outputContainer->Add(fhTimeENoCut);  
1421     
1422     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1423     fhTimeESPD->SetXTitle("E (GeV)");
1424     fhTimeESPD->SetYTitle("time (ns)");
1425     outputContainer->Add(fhTimeESPD);  
1426     
1427     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1428     fhTimeESPDMulti->SetXTitle("E (GeV)");
1429     fhTimeESPDMulti->SetYTitle("time (ns)");
1430     outputContainer->Add(fhTimeESPDMulti);  
1431     
1432     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50); 
1433     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1434     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1435     outputContainer->Add(fhTimeNPileUpVertSPD);  
1436     
1437     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 ); 
1438     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1439     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1440     outputContainer->Add(fhTimeNPileUpVertTrack);  
1441     
1442     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50); 
1443     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1444     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1445     outputContainer->Add(fhTimeNPileUpVertContributors);  
1446     
1447     fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50); 
1448     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1449     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1450     outputContainer->Add(fhTimePileUpMainVertexZDistance);  
1451     
1452     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50); 
1453     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1454     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1455     outputContainer->Add(fhTimePileUpMainVertexZDiamond);  
1456     
1457   }
1458   
1459   //Keep neutral meson selection histograms if requiered
1460   //Setting done in AliNeutralMesonSelection
1461   
1462   if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1463   {
1464     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
1465     
1466     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1467       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
1468     
1469     delete nmsHistos;
1470   }
1471   
1472   return outputContainer ;
1473   
1474 }
1475
1476 //_____________________________________________
1477 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
1478
1479   
1480   // Assign mc index depending on MC bit set
1481   
1482   if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
1483   {
1484     return kmcPi0 ;      
1485   }//pi0
1486   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
1487   {
1488     return kmcEta ; 
1489   }//eta          
1490   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
1491              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
1492   {
1493     return kmcConversion ; 
1494   }//conversion photon
1495   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
1496   {
1497     return kmcPhoton ; 
1498   }//photon   no conversion
1499   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
1500   {
1501     return kmcElectron ; 
1502   }//electron
1503   else 
1504   {
1505     return kmcHadron ; 
1506   }//other particles 
1507   
1508 }
1509
1510 //__________________________________________________________________
1511 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1, 
1512                                        AliAODPWG4Particle * photon2, 
1513                                        Int_t & label, Int_t & tag)
1514 {
1515   // Check the labels of pare in case mother was same pi0 or eta
1516   // Set the new AOD accordingly
1517   
1518   Int_t  label1 = photon1->GetLabel();
1519   Int_t  label2 = photon2->GetLabel();
1520   
1521   if(label1 < 0 || label2 < 0 ) return ;
1522   
1523   //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
1524   //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
1525   Int_t tag1 = photon1->GetTag();
1526   Int_t tag2 = photon2->GetTag();
1527   
1528   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1529   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && 
1530        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
1531      (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) && 
1532       GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
1533      )
1534   {
1535     
1536     //Check if pi0/eta mother is the same
1537     if(GetReader()->ReadStack())
1538     { 
1539       if(label1>=0)
1540       {
1541         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1542         label1 = mother1->GetFirstMother();
1543         //mother1 = GetMCStack()->Particle(label1);//pi0
1544       }
1545       if(label2>=0)
1546       {
1547         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1548         label2 = mother2->GetFirstMother();
1549         //mother2 = GetMCStack()->Particle(label2);//pi0
1550       }
1551     } // STACK
1552     else if(GetReader()->ReadAODMCParticles())
1553     {//&& (input > -1)){
1554       if(label1>=0)
1555       {
1556         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
1557         label1 = mother1->GetMother();
1558         //mother1 = GetMCStack()->Particle(label1);//pi0
1559       }
1560       if(label2>=0)
1561       {
1562         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
1563         label2 = mother2->GetMother();
1564         //mother2 = GetMCStack()->Particle(label2);//pi0
1565       }
1566     }// AOD
1567     
1568     //printf("mother1 %d, mother2 %d\n",label1,label2);
1569     if( label1 == label2 && label1>=0 )
1570     {
1571       
1572       label = label1;
1573       
1574       TLorentzVector mom1 = *(photon1->Momentum());
1575       TLorentzVector mom2 = *(photon2->Momentum());
1576       
1577       Double_t angle = mom2.Angle(mom1.Vect());
1578       Double_t mass  = (mom1+mom2).M();
1579       Double_t epair = (mom1+mom2).E();
1580       
1581       if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
1582       {
1583         fhMassPairMCPi0 ->Fill(epair,mass);
1584         fhAnglePairMCPi0->Fill(epair,angle);
1585         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1586       }
1587       else 
1588       {
1589         fhMassPairMCEta ->Fill(epair,mass);
1590         fhAnglePairMCEta->Fill(epair,angle);
1591         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
1592       }
1593       
1594     } // same label
1595   } // both from eta or pi0 decay
1596   
1597 }   
1598
1599 //____________________________________________________________________________
1600 void AliAnaPi0EbE::Init()
1601
1602   //Init
1603   //Do some checks
1604   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1605     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1606     abort();
1607   }
1608   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1609     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1610     abort();
1611   }
1612   
1613 }
1614
1615 //____________________________________________________________________________
1616 void AliAnaPi0EbE::InitParameters()
1617 {
1618   //Initialize the parameters of the analysis.  
1619   AddToHistogramsName("AnaPi0EbE_");
1620   
1621   fInputAODGammaConvName = "PhotonsCTS" ;
1622   fAnaType = kIMCalo ;
1623   fCalorimeter = "EMCAL" ;
1624   fMinDist  = 2.;
1625   fMinDist2 = 4.;
1626   fMinDist3 = 5.;
1627   
1628 }
1629
1630 //__________________________________________________________________
1631 void  AliAnaPi0EbE::MakeAnalysisFillAOD() 
1632 {
1633   //Do analysis and fill aods
1634   
1635   switch(fAnaType) 
1636   {
1637     case kIMCalo:
1638       MakeInvMassInCalorimeter();
1639       break;
1640       
1641     case kSSCalo:
1642       MakeShowerShapeIdentification();
1643       break;
1644       
1645     case kIMCaloTracks:
1646       MakeInvMassInCalorimeterAndCTS();
1647       break;
1648       
1649   }
1650 }
1651
1652 //____________________________________________
1653 void  AliAnaPi0EbE::MakeInvMassInCalorimeter() 
1654 {
1655   //Do analysis and fill aods
1656   //Search for the photon decay in calorimeters
1657   //Read photon list from AOD, produced in class AliAnaPhoton
1658   //Check if 2 photons have the mass of the pi0.
1659   
1660   TLorentzVector mom1;
1661   TLorentzVector mom2;
1662   TLorentzVector mom ;
1663
1664   Int_t tag   = 0;
1665   Int_t label = 0;
1666   
1667   if(!GetInputAODBranch()){
1668     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1669     abort();
1670   }
1671   
1672   //Get shower shape information of clusters
1673   TObjArray *clusters = 0;
1674   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1675   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
1676   
1677   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
1678     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1679     
1680     //Vertex cut in case of mixed events
1681     Int_t evtIndex1 = 0 ; 
1682     if(GetMixedEvent())
1683       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
1684     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
1685     mom1 = *(photon1->Momentum());
1686     
1687     //Get original cluster, to recover some information
1688     Int_t iclus = -1;
1689     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus); 
1690     
1691     if(!cluster1){
1692       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1693       return;
1694     }
1695     
1696     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1697     {
1698       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
1699       
1700       Int_t evtIndex2 = 0 ; 
1701       if(GetMixedEvent())
1702         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1703       
1704       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1705         continue ; 
1706       
1707       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
1708       
1709       mom2 = *(photon2->Momentum());
1710       
1711       //Get original cluster, to recover some information
1712       Int_t iclus2;
1713       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1); 
1714       
1715       if(!cluster2)
1716       {
1717         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1718         return;
1719       }
1720       
1721       Float_t e1    = photon1->E();      
1722       Float_t e2    = photon2->E();
1723       
1724       //Select clusters with good time window difference
1725       Float_t tof1  = cluster1->GetTOF()*1e9;;
1726       Float_t tof2  = cluster2->GetTOF()*1e9;;
1727       Double_t t12diff = tof1-tof2;
1728       fhEPairDiffTime->Fill(e1+e2,    t12diff);
1729       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1730       
1731       //Play with the MC stack if available
1732       if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1733
1734       // Check the invariant mass for different selection on the local maxima
1735       // Name of AOD method TO BE FIXED
1736       Int_t nMaxima1 = photon1->GetFiducialArea();
1737       Int_t nMaxima2 = photon2->GetFiducialArea();
1738       
1739       Double_t mass  = (mom1+mom2).M();
1740       Double_t epair = (mom1+mom2).E();
1741       
1742       if(nMaxima1==nMaxima2)
1743       {
1744         if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1745         else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1746         else                 fhMassPairLocMax[2]->Fill(epair,mass);
1747       }
1748       else if(nMaxima1==1 || nMaxima2==1)
1749       {
1750         if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
1751         else                             fhMassPairLocMax[4]->Fill(epair,mass); 
1752       }
1753       else  
1754         fhMassPairLocMax[5]->Fill(epair,mass);
1755       
1756       // combinations with SS axis cut and NLM cut
1757       if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass); 
1758       if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass); 
1759       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass); 
1760       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass); 
1761       
1762       //Skip events with too few or too many  NLM
1763       if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
1764       
1765       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
1766       
1767       //Mass of all pairs
1768       fhMass->Fill(epair,(mom1+mom2).M());
1769       
1770       //Select good pair (good phi, pt cuts, aperture and invariant mass)
1771       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1772       {
1773         if(GetDebug()>1) 
1774           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());
1775         
1776         //Fill some histograms about shower shape
1777         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1778         {
1779           FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1780           FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
1781         }
1782         
1783         // Tag both photons as decay
1784         photon1->SetTagged(kTRUE);
1785         photon2->SetTagged(kTRUE);
1786         
1787         fhPtDecay->Fill(photon1->Pt());
1788         fhEDecay ->Fill(photon1->E() );
1789         
1790         fhPtDecay->Fill(photon2->Pt());
1791         fhEDecay ->Fill(photon2->E() );
1792         
1793         //Create AOD for analysis
1794         mom = mom1+mom2;
1795                 
1796         //Mass of selected pairs
1797         fhSelectedMass->Fill(epair,mom.M());
1798         
1799         // Fill histograms to undertand pile-up before other cuts applied
1800         // Remember to relax time cuts in the reader
1801         FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);        
1802         
1803         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1804         
1805         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1806         pi0.SetDetector(photon1->GetDetector());
1807         
1808         // MC
1809         pi0.SetLabel(label);
1810         pi0.SetTag(tag);  
1811         
1812         //Set the indeces of the original caloclusters  
1813         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
1814         //pi0.SetInputFileIndex(input);
1815         
1816         AddAODParticle(pi0);
1817         
1818       }//pi0
1819       
1820     }//2n photon loop
1821     
1822   }//1st photon loop
1823   
1824   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");  
1825   
1826 }
1827
1828 //__________________________________________________
1829 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() 
1830 {
1831   //Do analysis and fill aods
1832   //Search for the photon decay in calorimeters
1833   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1834   //Check if 2 photons have the mass of the pi0.
1835   
1836   TLorentzVector mom1;
1837   TLorentzVector mom2;
1838   TLorentzVector mom ;
1839   Int_t tag   = 0;
1840   Int_t label = 0;
1841   Int_t evtIndex = 0;
1842   
1843   // Check calorimeter input
1844   if(!GetInputAODBranch()){
1845     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
1846     abort();
1847   }
1848   
1849   // Get the array with conversion photons
1850   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1851   if(!inputAODGammaConv) {
1852     
1853     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1854     
1855     if(!inputAODGammaConv) {
1856       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1857       
1858       return;
1859     }
1860   }  
1861   
1862   //Get shower shape information of clusters
1863   TObjArray *clusters = 0;
1864   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1865   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;  
1866   
1867   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
1868   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
1869   if(nCTS<=0 || nCalo <=0) 
1870   {
1871     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1872     return;
1873   }
1874   
1875   if(GetDebug() > 1)
1876     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1877   
1878   // Do the loop, first calo, second CTS
1879   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1880     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1881     mom1 = *(photon1->Momentum());
1882     
1883     //Get original cluster, to recover some information
1884     Int_t iclus = -1;
1885     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);     
1886     
1887     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1888       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
1889       if(GetMixedEvent())
1890         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1891       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
1892       
1893       mom2 = *(photon2->Momentum());
1894       
1895       Double_t mass  = (mom1+mom2).M();
1896       Double_t epair = (mom1+mom2).E();
1897       
1898       Int_t nMaxima = photon1->GetFiducialArea();
1899       if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1900       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1901       else                fhMassPairLocMax[2]->Fill(epair,mass);
1902       
1903       if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
1904       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
1905
1906       //Play with the MC stack if available
1907       if(IsDataMC())
1908       {
1909         Int_t   label2 = photon2->GetLabel();
1910         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
1911         
1912         HasPairSameMCMother(photon1, photon2, label, tag) ;
1913       }
1914       
1915       //Mass of selected pairs
1916       fhMass->Fill(epair,(mom1+mom2).M());
1917       
1918       //Select good pair (good phi, pt cuts, aperture and invariant mass)
1919       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1920       {
1921         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());
1922         
1923         //Fill some histograms about shower shape
1924         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1925         {
1926           FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1927         }        
1928         
1929         // Tag both photons as decay
1930         photon1->SetTagged(kTRUE);
1931         photon2->SetTagged(kTRUE);        
1932         
1933         fhPtDecay->Fill(photon1->Pt());
1934         fhEDecay ->Fill(photon1->E() );
1935         
1936         //Create AOD for analysis
1937         
1938         mom = mom1+mom2;
1939         
1940         //Mass of selected pairs
1941         fhSelectedMass->Fill(epair,mom.M());
1942         
1943         // Fill histograms to undertand pile-up before other cuts applied
1944         // Remember to relax time cuts in the reader
1945         if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
1946         
1947         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1948         
1949         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1950         pi0.SetDetector(photon1->GetDetector());
1951         
1952         // MC
1953         pi0.SetLabel(label);
1954         pi0.SetTag(tag);
1955         
1956         //Set the indeces of the original tracks or caloclusters  
1957         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1958         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
1959         //pi0.SetInputFileIndex(input);
1960         
1961         AddAODParticle(pi0);
1962         
1963       }//pi0
1964     }//2n photon loop
1965     
1966   }//1st photon loop
1967   
1968   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");  
1969   
1970 }
1971
1972
1973 //_________________________________________________
1974 void  AliAnaPi0EbE::MakeShowerShapeIdentification() 
1975 {
1976   //Search for pi0 in fCalorimeter with shower shape analysis 
1977   
1978   TObjArray * pl        = 0x0; 
1979   AliVCaloCells * cells = 0x0;
1980   //Select the Calorimeter of the photon
1981   if      (fCalorimeter == "PHOS" )
1982   {
1983     pl    = GetPHOSClusters();
1984     cells = GetPHOSCells();
1985   }
1986   else if (fCalorimeter == "EMCAL")
1987   {
1988     pl    = GetEMCALClusters();
1989     cells = GetEMCALCells();
1990   }
1991   
1992   if(!pl) 
1993   {
1994     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1995     return;
1996   }  
1997         
1998   TLorentzVector mom ;
1999   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2000   {
2001     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));        
2002     
2003     Int_t evtIndex = 0 ; 
2004     if (GetMixedEvent()) 
2005     {
2006       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
2007     }
2008     
2009     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
2010     
2011     //Get Momentum vector, 
2012     Double_t vertex[]={0,0,0};
2013     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2014     {
2015       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2016     }//Assume that come from vertex in straight line
2017     else
2018     {
2019       calo->GetMomentum(mom,vertex) ;
2020     }
2021           
2022     //If too small or big pt, skip it
2023     if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ; 
2024     
2025     //Check acceptance selection
2026     if(IsFiducialCutOn())
2027     {
2028       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2029       if(! in ) continue ;
2030     }
2031     
2032     if(GetDebug() > 1) 
2033       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());     
2034         
2035     //Check Distance to Bad channel, set bit.
2036     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2037     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2038     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
2039       continue ;
2040
2041     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2042
2043     //.......................................
2044     // TOF cut, BE CAREFUL WITH THIS CUT
2045     Double_t tof = calo->GetTOF()*1e9;
2046     if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
2047     
2048     //Play with the MC stack if available
2049     //Check origin of the candidates
2050     Int_t tag   = 0 ;
2051     if(IsDataMC())
2052     {
2053       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),0);
2054       //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2055       if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2056     }      
2057     
2058     //Skip matched clusters with tracks
2059     if(IsTrackMatched(calo, GetReader()->GetInputEvent())) 
2060     {
2061       FillRejectedClusterHistograms(mom,tag);
2062       continue ;
2063     }
2064         
2065     
2066     //Check PID
2067     //PID selection or bit setting
2068     Int_t    nMaxima = 0 ; 
2069     Double_t mass    = 0 , angle = 0;
2070     Double_t e1      = 0 , e2    = 0;
2071     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2072                                                                                    GetVertex(evtIndex),nMaxima,
2073                                                                                    mass,angle,e1,e2) ;   
2074     
2075     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
2076   
2077         
2078     //Skip events with too few or too many  NLM
2079     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) 
2080     {
2081       FillRejectedClusterHistograms(mom,tag);
2082       continue ;
2083     }
2084     
2085     if(GetDebug() > 1)
2086       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
2087     
2088     //mass of all clusters
2089     fhMass->Fill(mom.E(),mass);
2090
2091     // Asymmetry of all clusters
2092     Float_t asy =-10;      
2093     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
2094     fhAsymmetry->Fill(mom.E(),asy);
2095     
2096     if(IsDataMC())
2097     {
2098       Int_t mcIndex = GetMCIndex(tag);
2099       fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
2100     }  
2101     
2102     // If cluster does not pass pid, not pi0/eta, skip it.
2103     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0) 
2104     { 
2105       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
2106       FillRejectedClusterHistograms(mom,tag);
2107       continue ;
2108     }   
2109     
2110     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)     
2111     { 
2112       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
2113       FillRejectedClusterHistograms(mom,tag);
2114       continue ;
2115     }   
2116     
2117     if(GetDebug() > 1) 
2118       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
2119                               mom.Pt(), idPartType);
2120     
2121     fhSelectedAsymmetry->Fill(mom.E(),asy);
2122
2123     if( fUseSplitAsyCut &&  GetCaloPID()->IsInPi0SplitAsymmetryRange(mom.E(),asy,nMaxima) )
2124     {
2125       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Too large asymmetry\n");
2126       FillRejectedClusterHistograms(mom,tag);
2127       continue ;
2128     }
2129     
2130     //Mass of selected pairs
2131     fhSelectedMass     ->Fill(mom.E(),mass);
2132
2133     //-----------------------
2134     //Create AOD for analysis
2135     
2136     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
2137     aodpi0.SetLabel(calo->GetLabel());
2138     
2139     //Set the indeces of the original caloclusters  
2140     aodpi0.SetCaloLabel(calo->GetID(),-1);
2141     aodpi0.SetDetector(fCalorimeter);
2142
2143     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
2144     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ; 
2145     else                         aodpi0.SetDistToBad(0) ;
2146     
2147     // Check if cluster is pi0 via cluster splitting
2148     aodpi0.SetIdentifiedParticleType(idPartType); 
2149         
2150     // Add number of local maxima to AOD, method name in AOD to be FIXED
2151     aodpi0.SetFiducialArea(nMaxima);
2152     
2153     aodpi0.SetTag(tag);
2154     
2155     //Fill some histograms about shower shape
2156     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2157     {
2158       FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
2159     }  
2160     
2161     // Fill histograms to undertand pile-up before other cuts applied
2162     // Remember to relax time cuts in the reader
2163     FillPileUpHistograms(calo->E(),calo->GetTOF()*1e9);
2164     
2165     //Add AOD with pi0 object to aod branch
2166     AddAODParticle(aodpi0);
2167     
2168   }//loop
2169   
2170   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");  
2171   
2172 }
2173 //______________________________________________
2174 void  AliAnaPi0EbE::MakeAnalysisFillHistograms() 
2175 {
2176   //Do analysis and fill histograms
2177   
2178   if(!GetOutputAODBranch())
2179   {
2180     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
2181     abort();
2182   }
2183   //Loop on stored AOD pi0
2184   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2185   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2186   
2187   for(Int_t iaod = 0; iaod < naod ; iaod++)
2188   {
2189     
2190     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2191     Int_t pdg = pi0->GetIdentifiedParticleType();
2192           
2193     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;              
2194     
2195     //Fill pi0 histograms 
2196     Float_t ener  = pi0->E();
2197     Float_t pt    = pi0->Pt();
2198     Float_t phi   = pi0->Phi();
2199     if(phi < 0) phi+=TMath::TwoPi();
2200     Float_t eta = pi0->Eta();
2201     
2202     fhPt     ->Fill(pt);
2203     fhE      ->Fill(ener);
2204     
2205     fhEEta   ->Fill(ener,eta);
2206     fhEPhi   ->Fill(ener,phi);
2207     fhEtaPhi ->Fill(eta,phi);
2208
2209     if(fFillPileUpHistograms)
2210     {
2211       if(GetReader()->IsPileUpFromSPD())               fhPtPi0PileUp[0]->Fill(pt);
2212       if(GetReader()->IsPileUpFromEMCal())             fhPtPi0PileUp[1]->Fill(pt);
2213       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPi0PileUp[2]->Fill(pt);
2214       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPi0PileUp[3]->Fill(pt);
2215       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPi0PileUp[4]->Fill(pt);
2216       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPi0PileUp[5]->Fill(pt);
2217       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
2218     }
2219
2220     
2221     if(IsDataMC())
2222     {
2223       Int_t tag     = pi0->GetTag();
2224       Int_t mcIndex = GetMCIndex(tag);
2225
2226       fhMCE  [mcIndex] ->Fill(ener);
2227       fhMCPt [mcIndex] ->Fill(pt);
2228       fhMCPhi[mcIndex] ->Fill(pt,phi);
2229       fhMCEta[mcIndex] ->Fill(pt,eta);
2230       
2231       if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
2232       {
2233         Float_t efracMC = 0;
2234         Int_t label = pi0->GetLabel();
2235         
2236         Bool_t ok = kFALSE;
2237         TLorentzVector mom   = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok); 
2238         if(!ok) continue;
2239         
2240         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2241         {
2242           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok); 
2243           if(grandmom.E() > 0 && ok) 
2244           {
2245             efracMC =  grandmom.E()/ener;
2246             fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
2247           }
2248         }        
2249         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
2250         {
2251           fhMCPi0DecayPt->Fill(pt);
2252           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2253           if(grandmom.E() > 0 && ok) 
2254           {
2255             efracMC =  mom.E()/grandmom.E();
2256             fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
2257           }
2258         }
2259         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2260         {
2261           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok); 
2262           if(grandmom.E() > 0 && ok) 
2263           {
2264             efracMC =  grandmom.E()/ener;
2265             fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
2266           }
2267         }        
2268         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
2269         {
2270           fhMCEtaDecayPt->Fill(pt);
2271           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok); 
2272           if(grandmom.E() > 0 && ok) 
2273           {
2274             efracMC =  mom.E()/grandmom.E();
2275             fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
2276           }
2277         }
2278         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2279         {
2280           fhMCOtherDecayPt->Fill(pt);
2281         }
2282         
2283       }
2284       
2285     }//Histograms with MC
2286     
2287   }// aod loop
2288   
2289 }
2290
2291 //__________________________________________________________________
2292 void AliAnaPi0EbE::Print(const Option_t * opt) const
2293 {
2294   //Print some relevant parameters set for the analysis
2295   if(! opt)
2296     return;
2297   
2298   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2299   AliAnaCaloTrackCorrBaseClass::Print("");
2300   printf("Analysis Type = %d \n",  fAnaType) ;
2301   if(fAnaType == kSSCalo){     
2302     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
2303     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
2304     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2305     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); 
2306   } 
2307   printf("    \n") ;
2308   
2309
2310
2311