Add histogram with fraction of energy in cluster with large time cells
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPhoton.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 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 //
18 // Class for the photon identification.
19 // Clusters from calorimeters are identified as photons
20 // and kept in the AOD. Few histograms produced.
21 // Produces input for other analysis classes like AliAnaPi0, 
22 // AliAnaParticleHadronCorrelation ... 
23 //
24 // -- Author: Gustavo Conesa (LNF-INFN) 
25 //////////////////////////////////////////////////////////////////////////////
26   
27   
28 // --- ROOT system --- 
29 #include <TH2F.h>
30 #include <TH3D.h>
31 #include <TClonesArray.h>
32 #include <TObjString.h>
33 #include "TParticle.h"
34 #include "TDatabasePDG.h"
35
36 // --- Analysis system --- 
37 #include "AliAnaPhoton.h" 
38 #include "AliCaloTrackReader.h"
39 #include "AliStack.h"
40 #include "AliCaloPID.h"
41 #include "AliMCAnalysisUtils.h"
42 #include "AliFiducialCut.h"
43 #include "AliVCluster.h"
44 #include "AliAODMCParticle.h"
45 #include "AliMixedEvent.h"
46 #include "AliAODEvent.h"
47 #include "AliESDEvent.h"
48
49 // --- Detectors --- 
50 #include "AliPHOSGeoUtils.h"
51 #include "AliEMCALGeometry.h"
52
53 ClassImp(AliAnaPhoton)
54   
55 //____________________________
56 AliAnaPhoton::AliAnaPhoton() : 
57     AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""), 
58     fMinDist(0.),                 fMinDist2(0.),                fMinDist3(0.), 
59     fRejectTrackMatch(0),         fFillTMHisto(kFALSE),
60     fTimeCutMin(-10000),          fTimeCutMax(10000),         
61     fNCellsCut(0),                
62     fNLMCutMin(-1),               fNLMCutMax(10), 
63     fFillSSHistograms(kFALSE),    fFillOnlySimpleSSHisto(1),   
64     fNOriginHistograms(8),        fNPrimaryHistograms(4),
65     fFillPileUpHistograms(0),
66     // Histograms
67     fhNCellsE(0),                 fhCellsE(0),   // Control histograms            
68     fhMaxCellDiffClusterE(0),     fhTimeE(0),    // Control histograms
69     fhEPhoton(0),                 fhPtPhoton(0),  
70     fhPhiPhoton(0),               fhEtaPhoton(0), 
71     fhEtaPhiPhoton(0),            fhEtaPhi05Photon(0),
72
73     // Shower shape histograms
74     fhNLocMax(0),
75     fhDispE(0),                   fhLam0E(0),                   fhLam1E(0), 
76     fhDispETRD(0),                fhLam0ETRD(0),                fhLam1ETRD(0),
77     fhDispETM(0),                 fhLam0ETM(0),                 fhLam1ETM(0), 
78     fhDispETMTRD(0),              fhLam0ETMTRD(0),              fhLam1ETMTRD(0),
79
80     fhNCellsLam0LowE(0),          fhNCellsLam1LowE(0),          fhNCellsDispLowE(0),  
81     fhNCellsLam0HighE(0),         fhNCellsLam1HighE(0),         fhNCellsDispHighE(0),
82
83     fhEtaLam0LowE(0),             fhPhiLam0LowE(0), 
84     fhEtaLam0HighE(0),            fhPhiLam0HighE(0), 
85     fhLam0DispLowE(0),            fhLam0DispHighE(0), 
86     fhLam1Lam0LowE(0),            fhLam1Lam0HighE(0), 
87     fhDispLam1LowE(0),            fhDispLam1HighE(0),
88     fhDispEtaE(0),                fhDispPhiE(0),
89     fhSumEtaE(0),                 fhSumPhiE(0),                 fhSumEtaPhiE(0),
90     fhDispEtaPhiDiffE(0),         fhSphericityE(0),
91     fhDispSumEtaDiffE(0),         fhDispSumPhiDiffE(0),
92
93     // MC histograms
94     fhMCPhotonELambda0NoOverlap(0),       fhMCPhotonELambda0TwoOverlap(0),      fhMCPhotonELambda0NOverlap(0),
95     // Embedding
96     fhEmbeddedSignalFractionEnergy(0),     
97     fhEmbedPhotonELambda0FullSignal(0),   fhEmbedPhotonELambda0MostlySignal(0),  
98     fhEmbedPhotonELambda0MostlyBkg(0),    fhEmbedPhotonELambda0FullBkg(0),       
99     fhEmbedPi0ELambda0FullSignal(0),      fhEmbedPi0ELambda0MostlySignal(0),    
100     fhEmbedPi0ELambda0MostlyBkg(0),       fhEmbedPi0ELambda0FullBkg(0),
101     // PileUp
102     fhTimeENoCut(0),                      fhTimeESPD(0),        fhTimeESPDMulti(0),
103     fhTimeNPileUpVertSPD(0),              fhTimeNPileUpVertTrack(0),
104     fhTimeNPileUpVertContributors(0),
105     fhTimePileUpMainVertexZDistance(0),   fhTimePileUpMainVertexZDiamond(0),
106     fhClusterMultSPDPileUp(),             fhClusterMultNoPileUp()
107  {
108   //default ctor
109   
110   for(Int_t i = 0; i < 14; i++)
111   {
112     fhMCPt     [i] = 0;
113     fhMCE      [i] = 0;
114     fhMCPhi    [i] = 0;
115     fhMCEta    [i] = 0;
116     fhMCDeltaE [i] = 0;                
117     fhMCDeltaPt[i] = 0;
118     fhMC2E     [i] = 0;              
119     fhMC2Pt    [i] = 0;
120   }
121   
122   for(Int_t i = 0; i < 7; i++)
123   {
124     fhPtPrimMC [i] = 0;
125     fhEPrimMC  [i] = 0;
126     fhPhiPrimMC[i] = 0;
127     fhYPrimMC  [i] = 0;
128     
129     fhPtPrimMCAcc [i] = 0;
130     fhEPrimMCAcc  [i] = 0;
131     fhPhiPrimMCAcc[i] = 0;
132     fhYPrimMCAcc  [i] = 0;
133     
134     fhDispEtaDispPhi[i] = 0;
135     fhLambda0DispPhi[i] = 0;
136     fhLambda0DispEta[i] = 0;
137     
138     fhPtPileUp       [i] = 0;
139     fhPtChargedPileUp[i] = 0;
140     fhPtPhotonPileUp [i] = 0;
141
142     fhLambda0PileUp       [i] = 0;
143     fhLambda0ChargedPileUp[i] = 0;
144     
145     fhClusterEFracLongTimePileUp  [i] = 0;
146     
147     fhClusterTimeDiffPileUp       [i] = 0;
148     fhClusterTimeDiffChargedPileUp[i] = 0;
149     fhClusterTimeDiffPhotonPileUp [i] = 0;
150
151     for(Int_t j = 0; j < 6; j++)
152     {
153       fhMCDispEtaDispPhi[i][j] = 0;
154       fhMCLambda0DispEta[i][j] = 0;
155       fhMCLambda0DispPhi[i][j] = 0;
156     }
157   }  
158   
159   for(Int_t i = 0; i < 6; i++)
160   {
161     fhMCELambda0    [i]                  = 0;
162     fhMCELambda1    [i]                  = 0;
163     fhMCEDispersion [i]                  = 0;
164     fhMCNCellsE     [i]                  = 0; 
165     fhMCMaxCellDiffClusterE[i]           = 0; 
166     fhLambda0DispEta[i]                  = 0;
167     fhLambda0DispPhi[i]                  = 0;
168
169     fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
170     fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
171     fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
172     fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
173     fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
174     fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
175     
176     fhMCEDispEta       [i]               = 0;
177     fhMCEDispPhi       [i]               = 0;
178     fhMCESumEtaPhi     [i]               = 0;
179     fhMCEDispEtaPhiDiff[i]               = 0;
180     fhMCESphericity    [i]               = 0;
181   }
182   
183    for(Int_t i = 0; i < 5; i++) 
184    {
185      fhClusterCuts[i] = 0;
186    }
187   
188    // Track matching residuals
189    for(Int_t i = 0; i < 2; i++)
190    {
191      fhTrackMatchedDEta[i] = 0;                fhTrackMatchedDPhi[i] = 0;         fhTrackMatchedDEtaDPhi[i] = 0;   
192      fhTrackMatchedDEtaTRD[i] = 0;             fhTrackMatchedDPhiTRD[i] = 0;          
193      fhTrackMatchedDEtaMCOverlap[i] = 0;       fhTrackMatchedDPhiMCOverlap[i] = 0;
194      fhTrackMatchedDEtaMCNoOverlap[i] = 0;     fhTrackMatchedDPhiMCNoOverlap[i] = 0;
195      fhTrackMatchedDEtaMCConversion[i] = 0;    fhTrackMatchedDPhiMCConversion[i] = 0;
196      fhTrackMatchedMCParticle[i] = 0;          fhTrackMatchedMCParticle[i] = 0;   
197      fhdEdx[i] = 0;                            fhEOverP[i] = 0;
198      fhEOverPTRD[i] = 0;
199    }
200    
201    for(Int_t i = 0; i < 4; i++) 
202    {
203      fhClusterMultSPDPileUp[i] = 0;
204      fhClusterMultNoPileUp [i] = 0;
205    }
206    
207   //Initialize parameters
208   InitParameters();
209
210 }
211
212 //_____________________________________________________________________________________________________
213 Bool_t  AliAnaPhoton::ClusterSelected(AliVCluster* calo, const TLorentzVector mom, const Int_t nMaxima) 
214 {
215   //Select clusters if they pass different cuts
216   
217   Float_t ptcluster = mom.Pt();
218   Float_t ecluster  = mom.E();
219   Float_t l0cluster = calo->GetM02();
220   
221   if(GetDebug() > 2)
222     printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f\n",
223            GetReader()->GetEventNumber(),
224            ecluster,ptcluster, mom.Phi()*TMath::RadToDeg(),mom.Eta());
225   
226   fhClusterCuts[1]->Fill(ecluster);
227   
228   //.......................................
229   //If too small or big energy, skip it
230   if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) return kFALSE ; 
231   
232   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
233   
234   fhClusterCuts[2]->Fill(ecluster);
235
236   if(fFillPileUpHistograms)
237   {
238     
239     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
240     AliVCaloCells* cells = 0;
241     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
242     else                        cells = GetPHOSCells();
243     
244     Float_t maxCellFraction = 0.;
245     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
246
247     Double_t tmax  = cells->GetCellTime(absIdMax);
248     GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
249     tmax*=1.e9;
250     
251     Bool_t okPhoton = kFALSE;
252     if( GetCaloPID()->GetIdentifiedParticleType(calo)== AliCaloPID::kPhoton) okPhoton = kTRUE;
253
254     Bool_t matched = IsTrackMatched(calo,GetReader()->GetInputEvent());
255     Float_t clusterLongTimeE = 0;
256     Float_t clusterOKTimeE   = 0;
257     //Loop on cells inside cluster
258     for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
259     {
260       Int_t absId  = calo->GetCellsAbsId()[ipos];
261       //if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
262       if(cells->GetCellAmplitude(absIdMax) > 0.1)
263       {
264         Double_t time  = cells->GetCellTime(absId);
265         Float_t  amp   = cells->GetCellAmplitude(absId);
266         Int_t    bc    = GetReader()->GetInputEvent()->GetBunchCrossNumber();
267         GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,time,cells);
268         time*=1e9;
269         
270         Float_t diff = (tmax-time);
271       
272         if(GetReader()->IsInTimeWindow(time,amp)) clusterOKTimeE   += amp;
273         else                                      clusterLongTimeE += amp;
274         
275         if(GetReader()->IsPileUpFromSPD())
276         {
277           fhClusterTimeDiffPileUp[0]->Fill(ecluster, diff);
278           if(!matched)
279           {
280             fhClusterTimeDiffChargedPileUp[0]->Fill(ecluster, diff);
281             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[0]->Fill(ecluster, diff);
282           }
283         }
284         
285         if(GetReader()->IsPileUpFromEMCal())
286         {
287           fhClusterTimeDiffPileUp[1]->Fill(ecluster, diff);
288           if(!matched)
289           {
290             fhClusterTimeDiffChargedPileUp[1]->Fill(ecluster, diff);
291             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[1]->Fill(ecluster, diff);
292           }
293         }
294
295         if(GetReader()->IsPileUpFromSPDOrEMCal())
296         {
297           fhClusterTimeDiffPileUp[2]->Fill(ecluster, diff);
298           if(!matched)
299           {
300             fhClusterTimeDiffChargedPileUp[2]->Fill(ecluster, diff);
301             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[2]->Fill(ecluster, diff);
302           }
303         }
304         
305         if(GetReader()->IsPileUpFromSPDAndEMCal())
306         {
307           fhClusterTimeDiffPileUp[3]->Fill(ecluster, diff);
308           if(!matched)
309           {
310             fhClusterTimeDiffChargedPileUp[3]->Fill(ecluster, diff);
311             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[3]->Fill(ecluster, diff);
312           }
313         }
314         
315         if(GetReader()->IsPileUpFromSPDAndNotEMCal())
316         {
317           fhClusterTimeDiffPileUp[4]->Fill(ecluster, diff);
318           if(!matched)
319           {
320             fhClusterTimeDiffChargedPileUp[4]->Fill(ecluster, diff);
321             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[4]->Fill(ecluster, diff);
322           }
323         }
324         
325         if(GetReader()->IsPileUpFromEMCalAndNotSPD())
326         {
327           fhClusterTimeDiffPileUp[5]->Fill(ecluster, diff);
328           if(!matched)
329           {
330             fhClusterTimeDiffChargedPileUp[5]->Fill(ecluster, diff);
331             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[5]->Fill(ecluster, diff);
332           }
333         }
334
335         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
336         {
337           fhClusterTimeDiffPileUp[6]->Fill(ecluster, diff);
338           if(!matched)
339           {
340             fhClusterTimeDiffChargedPileUp[6]->Fill(ecluster, diff);
341             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[6]->Fill(ecluster, diff);
342           }
343         }
344       }// Not max
345     }//loop
346     
347     Float_t frac = 0;
348     if(clusterLongTimeE+clusterOKTimeE > 0.001)
349       frac = clusterLongTimeE/(clusterLongTimeE+clusterOKTimeE);
350     //printf("E long %f, E OK %f, Fraction large time %f, E %f\n",clusterLongTimeE,clusterOKTimeE,frac,ecluster);
351     
352     if(GetReader()->IsPileUpFromSPD())               {fhPtPileUp[0]->Fill(ptcluster); fhLambda0PileUp[0]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[0]->Fill(ecluster,frac);}
353     if(GetReader()->IsPileUpFromEMCal())             {fhPtPileUp[1]->Fill(ptcluster); fhLambda0PileUp[1]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[1]->Fill(ecluster,frac);}
354     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtPileUp[2]->Fill(ptcluster); fhLambda0PileUp[2]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[2]->Fill(ecluster,frac);}
355     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtPileUp[3]->Fill(ptcluster); fhLambda0PileUp[3]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[3]->Fill(ecluster,frac);}
356     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtPileUp[4]->Fill(ptcluster); fhLambda0PileUp[4]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[4]->Fill(ecluster,frac);}
357     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtPileUp[5]->Fill(ptcluster); fhLambda0PileUp[5]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[5]->Fill(ecluster,frac);}
358     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(ptcluster); fhLambda0PileUp[6]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[6]->Fill(ecluster,frac);}
359         
360   }
361   
362   //.......................................
363   // TOF cut, BE CAREFUL WITH THIS CUT
364   Double_t tof = calo->GetTOF()*1e9;
365   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
366   
367   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
368   
369   fhClusterCuts[3]->Fill(ecluster);
370
371   //.......................................
372   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
373   
374   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
375   
376   fhClusterCuts[4]->Fill(ecluster);
377
378   if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
379   if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
380
381   fhClusterCuts[5]->Fill(ecluster);
382
383   //.......................................
384   //Check acceptance selection
385   if(IsFiducialCutOn())
386   {
387     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
388     if(! in ) return kFALSE ;
389   }
390   
391   if(GetDebug() > 2) printf("Fiducial cut passed \n");
392   
393   fhClusterCuts[6]->Fill(ecluster);
394
395   //.......................................
396   //Skip matched clusters with tracks
397   
398   // Fill matching residual histograms before PID cuts
399   if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,0);
400   
401   if(fRejectTrackMatch)
402   {
403     if(IsTrackMatched(calo,GetReader()->GetInputEvent())) 
404     {
405       if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
406       return kFALSE ;
407     }
408     else  
409       if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
410   }// reject matched clusters
411   
412   fhClusterCuts[7]->Fill(ecluster);
413
414   if(fFillPileUpHistograms)
415   {
416     if(GetReader()->IsPileUpFromSPD())               {fhPtChargedPileUp[0]->Fill(ptcluster); fhLambda0ChargedPileUp[0]->Fill(ecluster,l0cluster); }
417     if(GetReader()->IsPileUpFromEMCal())             {fhPtChargedPileUp[1]->Fill(ptcluster); fhLambda0ChargedPileUp[1]->Fill(ecluster,l0cluster); }
418     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtChargedPileUp[2]->Fill(ptcluster); fhLambda0ChargedPileUp[2]->Fill(ecluster,l0cluster); }
419     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtChargedPileUp[3]->Fill(ptcluster); fhLambda0ChargedPileUp[3]->Fill(ecluster,l0cluster); }
420     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtChargedPileUp[4]->Fill(ptcluster); fhLambda0ChargedPileUp[4]->Fill(ecluster,l0cluster); }
421     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtChargedPileUp[5]->Fill(ptcluster); fhLambda0ChargedPileUp[5]->Fill(ecluster,l0cluster); }
422     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtChargedPileUp[6]->Fill(ptcluster); fhLambda0ChargedPileUp[6]->Fill(ecluster,l0cluster); }
423   }
424   
425   //.......................................
426   //Check Distance to Bad channel, set bit.
427   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
428   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
429   if(distBad < fMinDist) 
430   {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
431     return kFALSE ;
432   }
433   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
434   
435   fhClusterCuts[8]->Fill(ecluster);
436   
437   if(GetDebug() > 0) 
438     printf("AliAnaPhoton::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f\n",
439            GetReader()->GetEventNumber(), 
440            ecluster, ptcluster,mom.Phi()*TMath::RadToDeg(),mom.Eta());
441   
442   //All checks passed, cluster selected
443   return kTRUE;
444     
445 }
446
447 //___________________________________________
448 void AliAnaPhoton::FillAcceptanceHistograms()
449 {
450   //Fill acceptance histograms if MC data is available
451   
452   Double_t photonY   = -100 ;
453   Double_t photonE   = -1 ;
454   Double_t photonPt  = -1 ;
455   Double_t photonPhi =  100 ;
456   Double_t photonEta = -1 ;
457
458   Int_t    pdg       =  0 ;
459   Int_t    tag       =  0 ;
460   Int_t    mcIndex   =  0 ;
461   Bool_t   inacceptance = kFALSE;
462
463   if(GetReader()->ReadStack())
464   {     
465     AliStack * stack = GetMCStack();
466     if(stack)
467     {
468       for(Int_t i=0 ; i<stack->GetNtrack(); i++)
469       {
470         TParticle * prim = stack->Particle(i) ;
471         pdg = prim->GetPdgCode();
472         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
473         //                             prim->GetName(), prim->GetPdgCode());
474         
475         if(pdg == 22)
476         {
477           // Get tag of this particle photon from fragmentation, decay, prompt ...
478           tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
479           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
480           {
481             //A conversion photon from a hadron, skip this kind of photon
482             // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
483             // GetMCAnalysisUtils()->PrintMCTag(tag);
484             
485             return;
486           }
487           
488           //Get photon kinematics
489           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception          
490           
491           photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
492           photonE   = prim->Energy() ;
493           photonPt  = prim->Pt() ;
494           photonPhi = TMath::RadToDeg()*prim->Phi() ;
495           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
496           photonEta = prim->Eta() ;
497           
498           //Check if photons hit the Calorimeter
499           TLorentzVector lv;
500           prim->Momentum(lv);
501           inacceptance = kFALSE;
502           if     (fCalorimeter == "PHOS")
503           {
504             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
505             {
506               Int_t mod ;
507               Double_t x,z ;
508               if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x)) 
509                 inacceptance = kTRUE;
510               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
511             }
512             else
513             {
514               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
515                 inacceptance = kTRUE ;
516               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
517             }
518           }        
519           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
520           {
521             if(GetEMCALGeometry())
522             {
523               Int_t absID=0;
524               
525               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
526               
527               if( absID >= 0) 
528                 inacceptance = kTRUE;
529               
530               //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2)) 
531               //                    inacceptance = kTRUE;
532               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
533             }
534             else
535             {
536               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
537                 inacceptance = kTRUE ;
538               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
539             }
540           }       //In EMCAL
541           
542           //Fill histograms
543           fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
544           if(TMath::Abs(photonY) < 1.0)
545           {
546             fhEPrimMC  [kmcPPhoton]->Fill(photonE ) ;
547             fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
548             fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
549             fhYPrimMC  [kmcPPhoton]->Fill(photonE , photonEta) ;
550           }
551           if(inacceptance)
552           {
553             fhEPrimMCAcc  [kmcPPhoton]->Fill(photonE ) ;
554             fhPtPrimMCAcc [kmcPPhoton]->Fill(photonPt) ;
555             fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
556             fhYPrimMCAcc  [kmcPPhoton]->Fill(photonE , photonY) ;
557           }//Accepted
558           
559           //Origin of photon
560           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
561           {
562             mcIndex = kmcPPrompt;
563           }
564           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
565           {
566             mcIndex = kmcPFragmentation ;
567           }
568           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
569           {
570             mcIndex = kmcPISR;
571           }
572           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
573           {
574             mcIndex = kmcPPi0Decay;
575           }
576           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || 
577                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
578           {
579             mcIndex = kmcPOtherDecay;
580           }
581           else if(fhEPrimMC[kmcPOther])
582           {
583             mcIndex = kmcPOther;
584           }//Other origin
585           
586           fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
587           if(TMath::Abs(photonY) < 1.0)
588           {
589             fhEPrimMC  [mcIndex]->Fill(photonE ) ;
590             fhPtPrimMC [mcIndex]->Fill(photonPt) ;
591             fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
592             fhYPrimMC  [mcIndex]->Fill(photonE , photonEta) ;
593           }
594           if(inacceptance)
595           {
596             fhEPrimMCAcc  [mcIndex]->Fill(photonE ) ;
597             fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
598             fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
599             fhYPrimMCAcc  [mcIndex]->Fill(photonE , photonY) ;
600           }//Accepted
601           
602         }// Primary photon 
603       }//loop on primaries      
604     }//stack exists and data is MC
605   }//read stack
606   else if(GetReader()->ReadAODMCParticles())
607   {
608     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
609     if(mcparticles)
610     {
611       Int_t nprim = mcparticles->GetEntriesFast();
612       
613       for(Int_t i=0; i < nprim; i++)
614       {
615         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);   
616         
617         pdg = prim->GetPdgCode();
618         
619         if(pdg == 22)
620         {
621           // Get tag of this particle photon from fragmentation, decay, prompt ...
622           tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
623           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
624           {
625             //A conversion photon from a hadron, skip this kind of photon
626             //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
627             //            GetMCAnalysisUtils()->PrintMCTag(tag);
628             
629             return;
630           }
631           
632           //Get photon kinematics
633           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception      
634           
635           photonY   = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
636           photonE   = prim->E() ;
637           photonPt  = prim->Pt() ;
638           photonPhi = prim->Phi() ;
639           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
640           photonEta = prim->Eta() ;
641           
642           //Check if photons hit the Calorimeter
643           TLorentzVector lv;
644           lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
645           inacceptance = kFALSE;
646           if     (fCalorimeter == "PHOS")
647           {
648             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
649             {
650               Int_t mod ;
651               Double_t x,z ;
652               Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
653               if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
654                 inacceptance = kTRUE;
655               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
656             }
657             else
658             {
659               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
660                 inacceptance = kTRUE ;
661               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
662             }
663           }        
664           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
665           {
666             if(GetEMCALGeometry())
667             {
668               Int_t absID=0;
669               
670               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
671               
672               if( absID >= 0) 
673                 inacceptance = kTRUE;
674               
675               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
676             }
677             else
678             {
679               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
680                 inacceptance = kTRUE ;
681               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
682             }
683           }       //In EMCAL
684           
685           //Fill histograms
686           
687           fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
688           if(TMath::Abs(photonY) < 1.0)
689           {
690             fhEPrimMC  [kmcPPhoton]->Fill(photonE ) ;
691             fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
692             fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
693             fhYPrimMC[kmcPPhoton]->Fill(photonE , photonEta) ;
694           }
695           
696           if(inacceptance)
697           {
698             fhEPrimMCAcc[kmcPPhoton]  ->Fill(photonE ) ;
699             fhPtPrimMCAcc[kmcPPhoton] ->Fill(photonPt) ;
700             fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
701             fhYPrimMCAcc[kmcPPhoton]  ->Fill(photonE , photonY) ;
702           }//Accepted
703           
704           //Origin of photon
705           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
706           {
707             mcIndex = kmcPPrompt;
708           }
709           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
710           {
711             mcIndex = kmcPFragmentation ;
712           }
713           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
714           {
715             mcIndex = kmcPISR;
716           }
717           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
718           {
719             mcIndex = kmcPPi0Decay;
720           }
721           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || 
722                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
723           {
724             mcIndex = kmcPOtherDecay;
725           }
726           else if(fhEPrimMC[kmcPOther])
727           {
728             mcIndex = kmcPOther;
729           }//Other origin
730           
731           fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
732           if(TMath::Abs(photonY) < 1.0)
733           {
734             fhEPrimMC  [mcIndex]->Fill(photonE ) ;
735             fhPtPrimMC [mcIndex]->Fill(photonPt) ;
736             fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
737             fhYPrimMC  [mcIndex]->Fill(photonE , photonEta) ;
738           }
739           if(inacceptance)
740           {
741             fhEPrimMCAcc  [mcIndex]->Fill(photonE ) ;
742             fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
743             fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
744             fhYPrimMCAcc  [mcIndex]->Fill(photonE , photonY) ;
745           }//Accepted
746                     
747         }// Primary photon 
748       }//loop on primaries      
749       
750     }//kmc array exists and data is MC
751   }     // read AOD MC
752 }
753
754 //___________________________________________________________________
755 void AliAnaPhoton::FillPileUpHistogramsPerEvent(TObjArray * clusters) 
756 {
757   // Fill some histograms per event to understand pile-up
758   // Open the time cut in the reader to be more meaningful
759   
760   if(!fFillPileUpHistograms) return;
761     
762   // Loop on clusters, get the maximum energy cluster as reference
763   Int_t nclusters = clusters->GetEntriesFast();
764   Int_t   idMax = 0; 
765   Float_t  eMax = 0;
766   Float_t  tMax = 0;
767   for(Int_t iclus = 0; iclus < nclusters ; iclus++)
768   {
769           AliVCluster * clus =  (AliVCluster*) (clusters->At(iclus));   
770     if(clus->E() > eMax && TMath::Abs(clus->GetTOF()*1e9) < 20)
771     {
772       eMax  = clus->E();
773       tMax  = clus->GetTOF()*1e9;
774       idMax = iclus;
775     }
776   }
777
778   if(eMax < 5) return;
779   
780   // Loop again on clusters to compare this max cluster t and the rest of the clusters, if E > 0.3
781   Int_t n20  = 0;
782   Int_t n40  = 0;
783   Int_t n    = 0;
784   Int_t nOK  = 0;
785
786   for(Int_t iclus = 0; iclus < nclusters ; iclus++)
787   {
788           AliVCluster * clus =  (AliVCluster*) (clusters->At(iclus));   
789     
790     if(clus->E() < 0.3 || iclus==idMax) continue;
791     
792     Float_t tdiff = TMath::Abs(tMax-clus->GetTOF()*1e9);
793     n++;
794     if(tdiff < 20) nOK++;
795     else
796     {
797       n20++;
798       if(tdiff > 40 ) n40++;
799     }
800   }
801   
802   // Check pile-up and fill histograms depending on the different cluster multiplicities
803   if(GetReader()->IsPileUpFromSPD())
804   {    
805     fhClusterMultSPDPileUp[0]->Fill(eMax,n  );
806     fhClusterMultSPDPileUp[1]->Fill(eMax,nOK);
807     fhClusterMultSPDPileUp[2]->Fill(eMax,n20);
808     fhClusterMultSPDPileUp[3]->Fill(eMax,n40);
809   }
810   else 
811   {
812     fhClusterMultNoPileUp[0]->Fill(eMax,n  );
813     fhClusterMultNoPileUp[1]->Fill(eMax,nOK);
814     fhClusterMultNoPileUp[2]->Fill(eMax,n20);
815     fhClusterMultNoPileUp[3]->Fill(eMax,n40);    
816   }  
817   
818 }
819
820
821 //_________________________________________________________________________________________________
822 void AliAnaPhoton::FillPileUpHistograms(const Float_t energy, const Float_t pt, const Float_t time)
823 {
824   // Fill some histograms to understand pile-up
825   if(!fFillPileUpHistograms) return;
826   
827   //printf("E %f, time %f\n",energy,time);
828   AliVEvent * event = GetReader()->GetInputEvent();
829   
830   if(GetReader()->IsPileUpFromSPD())               fhPtPhotonPileUp[0]->Fill(pt);
831   if(GetReader()->IsPileUpFromEMCal())             fhPtPhotonPileUp[1]->Fill(pt);
832   if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPhotonPileUp[2]->Fill(pt);
833   if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPhotonPileUp[3]->Fill(pt);
834   if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPhotonPileUp[4]->Fill(pt);
835   if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPhotonPileUp[5]->Fill(pt);
836   if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPhotonPileUp[6]->Fill(pt);
837   
838   fhTimeENoCut->Fill(energy,time);
839   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
840   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
841   
842   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
843   
844   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
845   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
846   
847   // N pile up vertices
848   Int_t nVerticesSPD    = -1;
849   Int_t nVerticesTracks = -1;
850   
851   if      (esdEv)
852   {
853     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
854     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
855
856   }//ESD
857   else if (aodEv)
858   {
859     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
860     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
861   }//AOD
862   
863   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
864   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
865   
866   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n", 
867   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
868   
869   Int_t ncont = -1;
870   Float_t z1 = -1, z2 = -1;
871   Float_t diamZ = -1;
872   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
873   {
874     if      (esdEv)
875     {
876       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
877       ncont=pv->GetNContributors();
878       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
879       z2 = pv->GetZ();
880       diamZ = esdEv->GetDiamondZ();
881     }//ESD
882     else if (aodEv)
883     {
884       AliAODVertex *pv=aodEv->GetVertex(iVert);
885       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
886       ncont=pv->GetNContributors();
887       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
888       z2=pv->GetZ();
889       diamZ = aodEv->GetDiamondZ();
890     }// AOD
891     
892     Double_t distZ  = TMath::Abs(z2-z1);
893     diamZ  = TMath::Abs(z2-diamZ);
894
895     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
896     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
897     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
898
899   }// loop
900 }
901
902 //____________________________________________________________________________________
903 void  AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag)
904 {
905     //Fill cluster Shower Shape histograms
906   
907   if(!fFillSSHistograms || GetMixedEvent()) return;
908
909   Float_t energy  = cluster->E();
910   Int_t   ncells  = cluster->GetNCells();
911   Float_t lambda0 = cluster->GetM02();
912   Float_t lambda1 = cluster->GetM20();
913   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
914   
915   TLorentzVector mom;
916   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
917   {
918     cluster->GetMomentum(mom,GetVertex(0)) ;
919   }//Assume that come from vertex in straight line
920   else
921   {
922     Double_t vertex[]={0,0,0};
923     cluster->GetMomentum(mom,vertex) ;
924   }
925   
926   Float_t eta = mom.Eta();
927   Float_t phi = mom.Phi();
928   if(phi < 0) phi+=TMath::TwoPi();
929   
930   fhLam0E ->Fill(energy,lambda0);
931   fhLam1E ->Fill(energy,lambda1);
932   fhDispE ->Fill(energy,disp);
933     
934   if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
935   {
936     fhLam0ETRD->Fill(energy,lambda0);
937     fhLam1ETRD->Fill(energy,lambda1);
938     fhDispETRD->Fill(energy,disp);
939   }
940   
941   Float_t l0   = 0., l1   = 0.;
942   Float_t dispp= 0., dEta = 0., dPhi    = 0.; 
943   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
944   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
945   {
946     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
947                                                                                  l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
948     //printf("AliAnaPhoton::FillShowerShapeHistogram - l0 %2.6f, l1 %2.6f, disp %2.6f, dEta %2.6f, dPhi %2.6f, sEta %2.6f, sPhi %2.6f, sEtaPhi %2.6f \n",
949     //       l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
950     //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
951     //       disp, dPhi+dEta );
952     fhDispEtaE        -> Fill(energy,dEta);
953     fhDispPhiE        -> Fill(energy,dPhi);
954     fhSumEtaE         -> Fill(energy,sEta);
955     fhSumPhiE         -> Fill(energy,sPhi);
956     fhSumEtaPhiE      -> Fill(energy,sEtaPhi);
957     fhDispEtaPhiDiffE -> Fill(energy,dPhi-dEta);
958     if(dEta+dPhi>0)fhSphericityE     -> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
959     if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.));
960     if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.));  
961     
962     Int_t ebin = -1;
963     if      (energy < 2 ) ebin = 0;
964     else if (energy < 4 ) ebin = 1;
965     else if (energy < 6 ) ebin = 2;
966     else if (energy < 10) ebin = 3;
967     else if (energy < 15) ebin = 4;  
968     else if (energy < 20) ebin = 5;  
969     else                  ebin = 6;  
970     
971     fhDispEtaDispPhi[ebin]->Fill(dEta   ,dPhi);
972     fhLambda0DispEta[ebin]->Fill(lambda0,dEta);
973     fhLambda0DispPhi[ebin]->Fill(lambda0,dPhi);
974     
975   }
976   
977   // if track-matching was of, check effect of track-matching residual cut 
978   
979   if(!fRejectTrackMatch)
980   {
981     Float_t dZ  = cluster->GetTrackDz();
982     Float_t dR  = cluster->GetTrackDx();
983     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
984     {
985       dR = 2000., dZ = 2000.;
986       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
987     }   
988     
989     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
990     {
991       fhLam0ETM ->Fill(energy,lambda0);
992       fhLam1ETM ->Fill(energy,lambda1);
993       fhDispETM ->Fill(energy,disp);
994       
995       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
996       {
997         fhLam0ETMTRD->Fill(energy,lambda0);
998         fhLam1ETMTRD->Fill(energy,lambda1);
999         fhDispETMTRD->Fill(energy,disp);
1000       }
1001     }
1002   }// if track-matching was of, check effect of matching residual cut  
1003   
1004   
1005   if(!fFillOnlySimpleSSHisto){
1006     if(energy < 2)
1007     {
1008       fhNCellsLam0LowE ->Fill(ncells,lambda0);
1009       fhNCellsLam1LowE ->Fill(ncells,lambda1);
1010       fhNCellsDispLowE ->Fill(ncells,disp);
1011       
1012       fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
1013       fhLam0DispLowE  ->Fill(lambda0,disp);
1014       fhDispLam1LowE  ->Fill(disp,lambda1);
1015       fhEtaLam0LowE   ->Fill(eta,lambda0);
1016       fhPhiLam0LowE   ->Fill(phi,lambda0);  
1017     }
1018     else 
1019     {
1020       fhNCellsLam0HighE ->Fill(ncells,lambda0);
1021       fhNCellsLam1HighE ->Fill(ncells,lambda1);
1022       fhNCellsDispHighE ->Fill(ncells,disp);
1023       
1024       fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
1025       fhLam0DispHighE  ->Fill(lambda0,disp);
1026       fhDispLam1HighE  ->Fill(disp,lambda1);
1027       fhEtaLam0HighE   ->Fill(eta, lambda0);
1028       fhPhiLam0HighE   ->Fill(phi, lambda0);
1029     }
1030   }
1031   
1032   if(IsDataMC())
1033   {
1034     AliVCaloCells* cells = 0;
1035     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1036     else                        cells = GetPHOSCells();
1037     
1038     //Fill histograms to check shape of embedded clusters
1039     Float_t fraction = 0;
1040     if(GetReader()->IsEmbeddedClusterSelectionOn())
1041     {//Only working for EMCAL
1042       Float_t clusterE = 0; // recalculate in case corrections applied.
1043       Float_t cellE    = 0;
1044       for(Int_t icell  = 0; icell < cluster->GetNCells(); icell++)
1045       {
1046         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
1047         clusterE+=cellE;  
1048         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
1049       }
1050       
1051       //Fraction of total energy due to the embedded signal
1052       fraction/=clusterE;
1053       
1054       if(GetDebug() > 1 ) 
1055         printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
1056       
1057       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
1058       
1059     }  // embedded fraction    
1060     
1061     // Get the fraction of the cluster energy that carries the cell with highest energy
1062     Int_t   absID           =-1 ;
1063     Float_t maxCellFraction = 0.;
1064     
1065     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
1066     
1067     // Check the origin and fill histograms
1068     
1069     Int_t mcIndex = -1;
1070     
1071     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)     && 
1072        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
1073        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)        &&
1074        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1075     {
1076       mcIndex = kmcssPhoton ;
1077
1078       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1079       {
1080         //Check particle overlaps in cluster
1081         
1082         // Compare the primary depositing more energy with the rest, 
1083         // if no photon/electron as comon ancestor (conversions), count as other particle
1084         Int_t ancPDG = 0, ancStatus = -1;
1085         TLorentzVector momentum; TVector3 prodVertex;
1086         Int_t ancLabel = 0;
1087         Int_t noverlaps = 1;      
1088         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) 
1089         {
1090           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], 
1091                                                                GetReader(),ancPDG,ancStatus,momentum,prodVertex);
1092           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
1093         }
1094         //printf("N overlaps %d \n",noverlaps);
1095         
1096         if(noverlaps == 1)
1097         {
1098           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
1099         }
1100         else if(noverlaps == 2)
1101         {        
1102           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
1103         }
1104         else if(noverlaps > 2)
1105         {          
1106           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
1107         }
1108         else 
1109         {
1110           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
1111         }
1112       }//No embedding
1113       
1114       //Fill histograms to check shape of embedded clusters
1115       if(GetReader()->IsEmbeddedClusterSelectionOn())
1116       {
1117         if     (fraction > 0.9) 
1118         {
1119           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
1120         }
1121         else if(fraction > 0.5)
1122         {
1123           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
1124         }
1125         else if(fraction > 0.1)
1126         { 
1127           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
1128         }
1129         else
1130         {
1131           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
1132         }
1133       } // embedded
1134       
1135     }//photon   no conversion
1136     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
1137     {
1138       mcIndex = kmcssElectron ;
1139     }//electron
1140     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) && 
1141                GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
1142     {
1143       mcIndex = kmcssConversion ;
1144     }//conversion photon
1145     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  )
1146     {
1147       mcIndex = kmcssPi0 ;
1148       
1149       //Fill histograms to check shape of embedded clusters
1150       if(GetReader()->IsEmbeddedClusterSelectionOn())
1151       {
1152         if     (fraction > 0.9) 
1153         {
1154           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
1155         }
1156         else if(fraction > 0.5)
1157         {
1158           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
1159         }
1160         else if(fraction > 0.1)
1161         { 
1162           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
1163         }
1164         else
1165         {
1166           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
1167         }
1168       } // embedded      
1169       
1170     }//pi0
1171     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  )
1172     {
1173       mcIndex = kmcssEta ;
1174     }//eta    
1175     else 
1176     {
1177       mcIndex = kmcssOther ; 
1178     }//other particles 
1179     
1180     fhMCELambda0           [mcIndex]->Fill(energy, lambda0);
1181     fhMCELambda1           [mcIndex]->Fill(energy, lambda1);
1182     fhMCEDispersion        [mcIndex]->Fill(energy, disp);
1183     fhMCNCellsE            [mcIndex]->Fill(energy, ncells);
1184     fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction);
1185     
1186     if(!fFillOnlySimpleSSHisto)
1187     {
1188       if     (energy < 2.)
1189       {
1190         fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction);
1191         fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells,  maxCellFraction);
1192       }
1193       else if(energy < 6.)
1194       {
1195         fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction);
1196         fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells,  maxCellFraction);
1197       }
1198       else
1199       {
1200         fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction);
1201         fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells,  maxCellFraction);
1202       }
1203       
1204       if(fCalorimeter == "EMCAL")
1205       {
1206         fhMCEDispEta        [mcIndex]-> Fill(energy,dEta);
1207         fhMCEDispPhi        [mcIndex]-> Fill(energy,dPhi);
1208         fhMCESumEtaPhi      [mcIndex]-> Fill(energy,sEtaPhi);
1209         fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy,dPhi-dEta);
1210         if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));  
1211         
1212         Int_t ebin = -1;
1213         if      (energy < 2 ) ebin = 0;
1214         else if (energy < 4 ) ebin = 1;
1215         else if (energy < 6 ) ebin = 2;
1216         else if (energy < 10) ebin = 3;
1217         else if (energy < 15) ebin = 4;  
1218         else if (energy < 20) ebin = 5;  
1219         else                  ebin = 6;  
1220         
1221         fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta   ,dPhi);
1222         fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0,dEta);
1223         fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0,dPhi);      
1224       }
1225     }
1226   }//MC data
1227   
1228 }
1229
1230 //__________________________________________________________________________
1231 void AliAnaPhoton::FillTrackMatchingResidualHistograms(AliVCluster* cluster, 
1232                                                        const Int_t cut)
1233 {
1234   // If selected, fill histograms with residuals of matched clusters, help to define track matching cut
1235   // Residual filled for different cuts 0 (No cut), after 1 PID cut
1236     
1237   Float_t dZ = cluster->GetTrackDz();
1238   Float_t dR = cluster->GetTrackDx();
1239   
1240   if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1241   {
1242     dR = 2000., dZ = 2000.;
1243     GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1244   }   
1245     
1246   if(fhTrackMatchedDEta[cut] && TMath::Abs(dR) < 999)
1247   {
1248     fhTrackMatchedDEta[cut]->Fill(cluster->E(),dZ);
1249     fhTrackMatchedDPhi[cut]->Fill(cluster->E(),dR);
1250     
1251     if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dZ,dR);
1252     
1253     Int_t nSMod = GetModuleNumber(cluster);
1254     
1255     if(fCalorimeter=="EMCAL" &&  nSMod > 5)
1256     {
1257       fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(),dZ);
1258       fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(),dR);
1259     }
1260     
1261     // Check dEdx and E/p of matched clusters
1262     
1263     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1264     {      
1265       
1266       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1267       
1268       if(track) 
1269       {
1270         
1271         Float_t dEdx   = track->GetTPCsignal();
1272         Float_t eOverp = cluster->E()/track->P();
1273         
1274         fhdEdx[cut]  ->Fill(cluster->E(), dEdx);
1275         fhEOverP[cut]->Fill(cluster->E(), eOverp);
1276         
1277         if(fCalorimeter=="EMCAL" && nSMod > 5)
1278           fhEOverPTRD[cut]->Fill(cluster->E(), eOverp);
1279         
1280         
1281       }
1282       else
1283           printf("AliAnaPhoton::FillTrackMatchingResidualHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1284       
1285       
1286       
1287       if(IsDataMC())
1288       {
1289         
1290         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(),GetReader(), 0);
1291         
1292         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
1293         {
1294           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1295                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5 );
1296           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5 );
1297           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5 );
1298           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5 );
1299           
1300           // Check if several particles contributed to cluster and discard overlapped mesons
1301           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) || 
1302              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1303           {
1304             if(cluster->GetNLabels()==1)
1305             {
1306               fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(),dZ);
1307               fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(),dR);
1308             }
1309             else 
1310             {
1311               fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(),dZ);
1312               fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(),dR);
1313             }
1314             
1315           }// Check overlaps
1316           
1317         }
1318         else
1319         {
1320           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1321                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5 );
1322           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5 );
1323           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5 );
1324           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5 );
1325           
1326           // Check if several particles contributed to cluster
1327           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) || 
1328              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1329           {
1330             fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(),dZ);
1331             fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(),dR);
1332             
1333           }// Check overlaps          
1334           
1335         }
1336         
1337       } // MC 
1338       
1339     } // residuals window
1340     
1341   } // Small residual
1342   
1343 }
1344
1345 //___________________________________________
1346 TObjString *  AliAnaPhoton::GetAnalysisCuts()
1347 {       
1348   //Save parameters used for analysis
1349   TString parList ; //this will be list of parameters used for this analysis.
1350   const Int_t buffersize = 255;
1351   char onePar[buffersize] ;
1352   
1353   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
1354   parList+=onePar ;     
1355   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1356   parList+=onePar ;
1357   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
1358   parList+=onePar ;
1359   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
1360   parList+=onePar ;
1361   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
1362   parList+=onePar ;
1363   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
1364   parList+=onePar ;  
1365   
1366   //Get parameters set in base class.
1367   parList += GetBaseParametersList() ;
1368   
1369   //Get parameters set in PID class.
1370   parList += GetCaloPID()->GetPIDParametersList() ;
1371   
1372   //Get parameters set in FiducialCut class (not available yet)
1373   //parlist += GetFidCut()->GetFidCutParametersList() 
1374   
1375   return new TObjString(parList) ;
1376 }
1377
1378 //________________________________________________________________________
1379 TList *  AliAnaPhoton::GetCreateOutputObjects()
1380 {  
1381   // Create histograms to be saved in output file and 
1382   // store them in outputContainer
1383   TList * outputContainer = new TList() ; 
1384   outputContainer->SetName("PhotonHistos") ; 
1385         
1386   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin(); 
1387   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin(); 
1388   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();   
1389   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax   = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin   = GetHistogramRanges()->GetHistoShowerShapeMin();
1390   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax    = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin    = GetHistogramRanges()->GetHistoNClusterCellMin(); 
1391   Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();       
1392
1393   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();          
1394   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();          
1395   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1396   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();          
1397   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();          
1398   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1399   
1400   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         
1401   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();         
1402   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1403   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       
1404   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();       
1405   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1406   
1407   Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins for SS studies
1408   
1409   TString cut[] = {"Open","Reader","E","Time","NCells","NLM","Fidutial","Matching","Bad","PID"};
1410   for (Int_t i = 0; i < 10 ;  i++) 
1411   {
1412     fhClusterCuts[i] = new TH1F(Form("hCut_%d_%s", i, cut[i].Data()),
1413                                 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
1414                                 nptbins,ptmin,ptmax); 
1415     fhClusterCuts[i]->SetYTitle("dN/dE ");
1416     fhClusterCuts[i]->SetXTitle("E (GeV)");
1417     outputContainer->Add(fhClusterCuts[i]) ;   
1418   }
1419   
1420   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax); 
1421   fhNCellsE->SetXTitle("E (GeV)");
1422   fhNCellsE->SetYTitle("# of cells in cluster");
1423   outputContainer->Add(fhNCellsE);    
1424   
1425   fhCellsE  = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax); 
1426   fhCellsE->SetXTitle("E_{cluster} (GeV)");
1427   fhCellsE->SetYTitle("E_{cell} (GeV)");
1428   outputContainer->Add(fhCellsE);    
1429   
1430   fhTimeE  = new TH2F ("hTimeE","time of cluster vs E of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1431   fhTimeE->SetXTitle("E (GeV)");
1432   fhTimeE->SetYTitle("time (ns)");
1433   outputContainer->Add(fhTimeE);  
1434   
1435   fhMaxCellDiffClusterE  = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1436                                     nptbins,ptmin,ptmax, 500,0,1.); 
1437   fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
1438   fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1439   outputContainer->Add(fhMaxCellDiffClusterE);  
1440   
1441   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax); 
1442   fhEPhoton->SetYTitle("N");
1443   fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
1444   outputContainer->Add(fhEPhoton) ;   
1445   
1446   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax); 
1447   fhPtPhoton->SetYTitle("N");
1448   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
1449   outputContainer->Add(fhPtPhoton) ; 
1450   
1451   fhPhiPhoton  = new TH2F
1452     ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
1453   fhPhiPhoton->SetYTitle("#phi (rad)");
1454   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1455   outputContainer->Add(fhPhiPhoton) ; 
1456   
1457   fhEtaPhoton  = new TH2F
1458     ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
1459   fhEtaPhoton->SetYTitle("#eta");
1460   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1461   outputContainer->Add(fhEtaPhoton) ;
1462   
1463   fhEtaPhiPhoton  = new TH2F
1464   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
1465   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
1466   fhEtaPhiPhoton->SetXTitle("#eta");
1467   outputContainer->Add(fhEtaPhiPhoton) ;
1468   if(GetMinPt() < 0.5)
1469   {
1470     fhEtaPhi05Photon  = new TH2F
1471     ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax); 
1472     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
1473     fhEtaPhi05Photon->SetXTitle("#eta");
1474     outputContainer->Add(fhEtaPhi05Photon) ;
1475   }
1476   
1477   fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
1478                        nptbins,ptmin,ptmax,10,0,10); 
1479   fhNLocMax ->SetYTitle("N maxima");
1480   fhNLocMax ->SetXTitle("E (GeV)");
1481   outputContainer->Add(fhNLocMax) ;  
1482   
1483   //Shower shape
1484   if(fFillSSHistograms)
1485   {
1486     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1487     fhLam0E->SetYTitle("#lambda_{0}^{2}");
1488     fhLam0E->SetXTitle("E (GeV)");
1489     outputContainer->Add(fhLam0E);  
1490     
1491     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1492     fhLam1E->SetYTitle("#lambda_{1}^{2}");
1493     fhLam1E->SetXTitle("E (GeV)");
1494     outputContainer->Add(fhLam1E);  
1495     
1496     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1497     fhDispE->SetYTitle("D^{2}");
1498     fhDispE->SetXTitle("E (GeV) ");
1499     outputContainer->Add(fhDispE);
1500
1501     if(!fRejectTrackMatch)
1502     {
1503       fhLam0ETM  = new TH2F ("hLam0ETM","#lambda_{0}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05,  |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1504       fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
1505       fhLam0ETM->SetXTitle("E (GeV)");
1506       outputContainer->Add(fhLam0ETM);  
1507       
1508       fhLam1ETM  = new TH2F ("hLam1ETM","#lambda_{1}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05,  |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1509       fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
1510       fhLam1ETM->SetXTitle("E (GeV)");
1511       outputContainer->Add(fhLam1ETM);  
1512       
1513       fhDispETM  = new TH2F ("hDispETM"," dispersion^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05,  |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1514       fhDispETM->SetYTitle("D^{2}");
1515       fhDispETM->SetXTitle("E (GeV) ");
1516       outputContainer->Add(fhDispETM);
1517     }
1518     
1519     if(fCalorimeter == "EMCAL")
1520     {
1521       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1522       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1523       fhLam0ETRD->SetXTitle("E (GeV)");
1524       outputContainer->Add(fhLam0ETRD);  
1525       
1526       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1527       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1528       fhLam1ETRD->SetXTitle("E (GeV)");
1529       outputContainer->Add(fhLam1ETRD);  
1530       
1531       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1532       fhDispETRD->SetYTitle("Dispersion^{2}");
1533       fhDispETRD->SetXTitle("E (GeV) ");
1534       outputContainer->Add(fhDispETRD);
1535       
1536       if(!fRejectTrackMatch)
1537       {
1538         fhLam0ETMTRD  = new TH2F ("hLam0ETMTRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05,  |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1539         fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
1540         fhLam0ETMTRD->SetXTitle("E (GeV)");
1541         outputContainer->Add(fhLam0ETMTRD);  
1542         
1543         fhLam1ETMTRD  = new TH2F ("hLam1ETMTRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05,  |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1544         fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
1545         fhLam1ETMTRD->SetXTitle("E (GeV)");
1546         outputContainer->Add(fhLam1ETMTRD);  
1547         
1548         fhDispETMTRD  = new TH2F ("hDispETMTRD"," dispersion^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05,  |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1549         fhDispETMTRD->SetYTitle("Dispersion^{2}");
1550         fhDispETMTRD->SetXTitle("E (GeV) ");
1551         outputContainer->Add(fhDispETMTRD);         
1552       } 
1553     } 
1554     
1555     if(!fFillOnlySimpleSSHisto)
1556     {
1557       fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1558       fhNCellsLam0LowE->SetXTitle("N_{cells}");
1559       fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1560       outputContainer->Add(fhNCellsLam0LowE);  
1561       
1562       fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1563       fhNCellsLam0HighE->SetXTitle("N_{cells}");
1564       fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1565       outputContainer->Add(fhNCellsLam0HighE);  
1566       
1567       fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1568       fhNCellsLam1LowE->SetXTitle("N_{cells}");
1569       fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1570       outputContainer->Add(fhNCellsLam1LowE);  
1571       
1572       fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1573       fhNCellsLam1HighE->SetXTitle("N_{cells}");
1574       fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1575       outputContainer->Add(fhNCellsLam1HighE);  
1576       
1577       fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1578       fhNCellsDispLowE->SetXTitle("N_{cells}");
1579       fhNCellsDispLowE->SetYTitle("D^{2}");
1580       outputContainer->Add(fhNCellsDispLowE);  
1581       
1582       fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1583       fhNCellsDispHighE->SetXTitle("N_{cells}");
1584       fhNCellsDispHighE->SetYTitle("D^{2}");
1585       outputContainer->Add(fhNCellsDispHighE);  
1586       
1587       fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax); 
1588       fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1589       fhEtaLam0LowE->SetXTitle("#eta");
1590       outputContainer->Add(fhEtaLam0LowE);  
1591       
1592       fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
1593       fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1594       fhPhiLam0LowE->SetXTitle("#phi");
1595       outputContainer->Add(fhPhiLam0LowE);  
1596       
1597       fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax); 
1598       fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1599       fhEtaLam0HighE->SetXTitle("#eta");
1600       outputContainer->Add(fhEtaLam0HighE);  
1601       
1602       fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
1603       fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1604       fhPhiLam0HighE->SetXTitle("#phi");
1605       outputContainer->Add(fhPhiLam0HighE);  
1606       
1607       fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1608       fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1609       fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1610       outputContainer->Add(fhLam1Lam0LowE);  
1611       
1612       fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1613       fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1614       fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1615       outputContainer->Add(fhLam1Lam0HighE);  
1616       
1617       fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1618       fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1619       fhLam0DispLowE->SetYTitle("D^{2}");
1620       outputContainer->Add(fhLam0DispLowE);  
1621       
1622       fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1623       fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1624       fhLam0DispHighE->SetYTitle("D^{2}");
1625       outputContainer->Add(fhLam0DispHighE);  
1626       
1627       fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1628       fhDispLam1LowE->SetXTitle("D^{2}");
1629       fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1630       outputContainer->Add(fhDispLam1LowE);  
1631       
1632       fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1633       fhDispLam1HighE->SetXTitle("D^{2}");
1634       fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1635       outputContainer->Add(fhDispLam1HighE);  
1636       
1637       if(fCalorimeter == "EMCAL")
1638       {
1639         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); 
1640         fhDispEtaE->SetXTitle("E (GeV)");
1641         fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
1642         outputContainer->Add(fhDispEtaE);     
1643         
1644         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); 
1645         fhDispPhiE->SetXTitle("E (GeV)");
1646         fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
1647         outputContainer->Add(fhDispPhiE);  
1648         
1649         fhSumEtaE  = new TH2F ("hSumEtaE","#delta^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
1650         fhSumEtaE->SetXTitle("E (GeV)");
1651         fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
1652         outputContainer->Add(fhSumEtaE);     
1653         
1654         fhSumPhiE  = new TH2F ("hSumPhiE","#delta^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",  
1655                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
1656         fhSumPhiE->SetXTitle("E (GeV)");
1657         fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
1658         outputContainer->Add(fhSumPhiE);  
1659         
1660         fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",  
1661                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
1662         fhSumEtaPhiE->SetXTitle("E (GeV)");
1663         fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
1664         outputContainer->Add(fhSumEtaPhiE);
1665         
1666         fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E", 
1667                                        nptbins,ptmin,ptmax,200, -10,10); 
1668         fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
1669         fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1670         outputContainer->Add(fhDispEtaPhiDiffE);    
1671         
1672         fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",  
1673                                    nptbins,ptmin,ptmax, 200, -1,1); 
1674         fhSphericityE->SetXTitle("E (GeV)");
1675         fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1676         outputContainer->Add(fhSphericityE);
1677         
1678         fhDispSumEtaDiffE  = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average vs E",  nptbins,ptmin,ptmax, 200,-0.01,0.01); 
1679         fhDispSumEtaDiffE->SetXTitle("E (GeV)");
1680         fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average");
1681         outputContainer->Add(fhDispSumEtaDiffE);     
1682         
1683         fhDispSumPhiDiffE  = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average vs E",  nptbins,ptmin,ptmax, 200,-0.01,0.01); 
1684         fhDispSumPhiDiffE->SetXTitle("E (GeV)");
1685         fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average");
1686         outputContainer->Add(fhDispSumPhiDiffE);  
1687         
1688         for(Int_t i = 0; i < 7; i++)
1689         {
1690           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]), 
1691                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
1692           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1693           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1694           outputContainer->Add(fhDispEtaDispPhi[i]); 
1695           
1696           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]), 
1697                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
1698           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1699           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1700           outputContainer->Add(fhLambda0DispEta[i]);       
1701           
1702           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]), 
1703                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
1704           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1705           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1706           outputContainer->Add(fhLambda0DispPhi[i]);         
1707         }
1708       }
1709     }
1710   } // Shower shape
1711   
1712   // Track Matching
1713   
1714   if(fFillTMHisto)
1715   {
1716     fhTrackMatchedDEta[0]  = new TH2F
1717     ("hTrackMatchedDEtaNoCut",
1718      "d#eta of cluster-track vs cluster energy, no photon cuts",
1719      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1720     fhTrackMatchedDEta[0]->SetYTitle("d#eta");
1721     fhTrackMatchedDEta[0]->SetXTitle("E_{cluster} (GeV)");
1722     
1723     fhTrackMatchedDPhi[0]  = new TH2F
1724     ("hTrackMatchedDPhiNoCut",
1725      "d#phi of cluster-track vs cluster energy, no photon cuts",
1726      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1727     fhTrackMatchedDPhi[0]->SetYTitle("d#phi (rad)");
1728     fhTrackMatchedDPhi[0]->SetXTitle("E_{cluster} (GeV)");
1729     
1730     fhTrackMatchedDEtaDPhi[0]  = new TH2F
1731     ("hTrackMatchedDEtaDPhiNoCut",
1732      "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1733      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax); 
1734     fhTrackMatchedDEtaDPhi[0]->SetYTitle("d#phi (rad)");
1735     fhTrackMatchedDEtaDPhi[0]->SetXTitle("d#eta");   
1736         
1737     fhdEdx[0]  = new TH2F ("hdEdxNoCut","matched track <dE/dx> vs cluster E, no photon cuts ", 
1738                            nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax); 
1739     fhdEdx[0]->SetXTitle("E (GeV)");
1740     fhdEdx[0]->SetYTitle("<dE/dx>");
1741     
1742     fhEOverP[0]  = new TH2F ("hEOverPNoCut","matched track E/p vs cluster E, no photon cuts ", 
1743                              nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax); 
1744     fhEOverP[0]->SetXTitle("E (GeV)");
1745     fhEOverP[0]->SetYTitle("E/p");
1746     
1747     outputContainer->Add(fhTrackMatchedDEta[0]) ; 
1748     outputContainer->Add(fhTrackMatchedDPhi[0]) ;
1749     outputContainer->Add(fhTrackMatchedDEtaDPhi[0]) ;
1750     outputContainer->Add(fhdEdx[0]);  
1751     outputContainer->Add(fhEOverP[0]);  
1752
1753     fhTrackMatchedDEta[1]  = new TH2F
1754     ("hTrackMatchedDEta",
1755      "d#eta of cluster-track vs cluster energy, no photon cuts",
1756      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1757     fhTrackMatchedDEta[1]->SetYTitle("d#eta");
1758     fhTrackMatchedDEta[1]->SetXTitle("E_{cluster} (GeV)");
1759     
1760     fhTrackMatchedDPhi[1]  = new TH2F
1761     ("hTrackMatchedDPhi",
1762      "d#phi of cluster-track vs cluster energy, no photon cuts",
1763      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1764     fhTrackMatchedDPhi[1]->SetYTitle("d#phi (rad)");
1765     fhTrackMatchedDPhi[1]->SetXTitle("E_{cluster} (GeV)");
1766     
1767     fhTrackMatchedDEtaDPhi[1]  = new TH2F
1768     ("hTrackMatchedDEtaDPhi",
1769      "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1770      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax); 
1771     fhTrackMatchedDEtaDPhi[1]->SetYTitle("d#phi (rad)");
1772     fhTrackMatchedDEtaDPhi[1]->SetXTitle("d#eta");   
1773     
1774     fhdEdx[1]  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", 
1775                            nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax); 
1776     fhdEdx[1]->SetXTitle("E (GeV)");
1777     fhdEdx[1]->SetYTitle("<dE/dx>");
1778     
1779     fhEOverP[1]  = new TH2F ("hEOverP","matched track E/p vs cluster E ", 
1780                              nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax); 
1781     fhEOverP[1]->SetXTitle("E (GeV)");
1782     fhEOverP[1]->SetYTitle("E/p");
1783     
1784     outputContainer->Add(fhTrackMatchedDEta[1]) ; 
1785     outputContainer->Add(fhTrackMatchedDPhi[1]) ;
1786     outputContainer->Add(fhTrackMatchedDEtaDPhi[1]) ;
1787     outputContainer->Add(fhdEdx[1]);  
1788     outputContainer->Add(fhEOverP[1]);      
1789     
1790     if(fCalorimeter=="EMCAL")
1791     {
1792       fhTrackMatchedDEtaTRD[0]  = new TH2F
1793       ("hTrackMatchedDEtaTRDNoCut",
1794        "d#eta of cluster-track vs cluster energy, SM behind TRD, no photon cuts",
1795        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1796       fhTrackMatchedDEtaTRD[0]->SetYTitle("d#eta");
1797       fhTrackMatchedDEtaTRD[0]->SetXTitle("E_{cluster} (GeV)");
1798       
1799       fhTrackMatchedDPhiTRD[0]  = new TH2F
1800       ("hTrackMatchedDPhiTRDNoCut",
1801        "d#phi of cluster-track vs cluster energy, SM behing TRD, no photon cuts",
1802        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1803       fhTrackMatchedDPhiTRD[0]->SetYTitle("d#phi (rad)");
1804       fhTrackMatchedDPhiTRD[0]->SetXTitle("E_{cluster} (GeV)");
1805             
1806       fhEOverPTRD[0]  = new TH2F ("hEOverPTRDNoCut","matched track E/p vs cluster E, behind TRD, no photon cuts ", 
1807                                   nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax); 
1808       fhEOverPTRD[0]->SetXTitle("E (GeV)");
1809       fhEOverPTRD[0]->SetYTitle("E/p");
1810       
1811       outputContainer->Add(fhTrackMatchedDEtaTRD[0]) ; 
1812       outputContainer->Add(fhTrackMatchedDPhiTRD[0]) ;
1813       outputContainer->Add(fhEOverPTRD[0]);  
1814       
1815       fhTrackMatchedDEtaTRD[1]  = new TH2F
1816       ("hTrackMatchedDEtaTRD",
1817        "d#eta of cluster-track vs cluster energy, SM behind TRD",
1818        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1819       fhTrackMatchedDEtaTRD[1]->SetYTitle("d#eta");
1820       fhTrackMatchedDEtaTRD[1]->SetXTitle("E_{cluster} (GeV)");
1821       
1822       fhTrackMatchedDPhiTRD[1]  = new TH2F
1823       ("hTrackMatchedDPhiTRD",
1824        "d#phi of cluster-track vs cluster energy, SM behing TRD",
1825        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1826       fhTrackMatchedDPhiTRD[1]->SetYTitle("d#phi (rad)");
1827       fhTrackMatchedDPhiTRD[1]->SetXTitle("E_{cluster} (GeV)");
1828       
1829       fhEOverPTRD[1]  = new TH2F ("hEOverPTRD","matched track E/p vs cluster E, behind TRD ", 
1830                                   nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax); 
1831       fhEOverPTRD[1]->SetXTitle("E (GeV)");
1832       fhEOverPTRD[1]->SetYTitle("E/p");
1833       
1834       outputContainer->Add(fhTrackMatchedDEtaTRD[1]) ; 
1835       outputContainer->Add(fhTrackMatchedDPhiTRD[1]) ;
1836       outputContainer->Add(fhEOverPTRD[1]);  
1837       
1838     }
1839     
1840     if(IsDataMC())
1841     {
1842       fhTrackMatchedDEtaMCNoOverlap[0]  = new TH2F
1843       ("hTrackMatchedDEtaMCNoOverlapNoCut",
1844        "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1845        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1846       fhTrackMatchedDEtaMCNoOverlap[0]->SetYTitle("d#eta");
1847       fhTrackMatchedDEtaMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1848       
1849       fhTrackMatchedDPhiMCNoOverlap[0]  = new TH2F
1850       ("hTrackMatchedDPhiMCNoOverlapNoCut",
1851        "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1852        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1853       fhTrackMatchedDPhiMCNoOverlap[0]->SetYTitle("d#phi (rad)");
1854       fhTrackMatchedDPhiMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1855       
1856       outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[0]) ; 
1857       outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[0]) ;
1858       
1859       fhTrackMatchedDEtaMCNoOverlap[1]  = new TH2F
1860       ("hTrackMatchedDEtaMCNoOverlap",
1861        "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1862        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1863       fhTrackMatchedDEtaMCNoOverlap[1]->SetYTitle("d#eta");
1864       fhTrackMatchedDEtaMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1865       
1866       fhTrackMatchedDPhiMCNoOverlap[1]  = new TH2F
1867       ("hTrackMatchedDPhiMCNoOverlap",
1868        "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1869        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1870       fhTrackMatchedDPhiMCNoOverlap[1]->SetYTitle("d#phi (rad)");
1871       fhTrackMatchedDPhiMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1872       
1873       outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[1]) ; 
1874       outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[1]) ;
1875       
1876       fhTrackMatchedDEtaMCOverlap[0]  = new TH2F
1877       ("hTrackMatchedDEtaMCOverlapNoCut",
1878        "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1879        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1880       fhTrackMatchedDEtaMCOverlap[0]->SetYTitle("d#eta");
1881       fhTrackMatchedDEtaMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1882       
1883       fhTrackMatchedDPhiMCOverlap[0]  = new TH2F
1884       ("hTrackMatchedDPhiMCOverlapNoCut",
1885        "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1886        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1887       fhTrackMatchedDPhiMCOverlap[0]->SetYTitle("d#phi (rad)");
1888       fhTrackMatchedDPhiMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1889       
1890       outputContainer->Add(fhTrackMatchedDEtaMCOverlap[0]) ; 
1891       outputContainer->Add(fhTrackMatchedDPhiMCOverlap[0]) ;
1892       
1893       fhTrackMatchedDEtaMCOverlap[1]  = new TH2F
1894       ("hTrackMatchedDEtaMCOverlap",
1895        "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1896        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1897       fhTrackMatchedDEtaMCOverlap[1]->SetYTitle("d#eta");
1898       fhTrackMatchedDEtaMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1899       
1900       fhTrackMatchedDPhiMCOverlap[1]  = new TH2F
1901       ("hTrackMatchedDPhiMCOverlap",
1902        "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1903        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1904       fhTrackMatchedDPhiMCOverlap[1]->SetYTitle("d#phi (rad)");
1905       fhTrackMatchedDPhiMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1906       
1907       outputContainer->Add(fhTrackMatchedDEtaMCOverlap[1]) ; 
1908       outputContainer->Add(fhTrackMatchedDPhiMCOverlap[1]) ;      
1909       
1910       fhTrackMatchedDEtaMCConversion[0]  = new TH2F
1911       ("hTrackMatchedDEtaMCConversionNoCut",
1912        "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1913        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1914       fhTrackMatchedDEtaMCConversion[0]->SetYTitle("d#eta");
1915       fhTrackMatchedDEtaMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
1916       
1917       fhTrackMatchedDPhiMCConversion[0]  = new TH2F
1918       ("hTrackMatchedDPhiMCConversionNoCut",
1919        "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1920        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1921       fhTrackMatchedDPhiMCConversion[0]->SetYTitle("d#phi (rad)");
1922       fhTrackMatchedDPhiMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
1923       
1924       outputContainer->Add(fhTrackMatchedDEtaMCConversion[0]) ; 
1925       outputContainer->Add(fhTrackMatchedDPhiMCConversion[0]) ;
1926        
1927       
1928       fhTrackMatchedDEtaMCConversion[1]  = new TH2F
1929       ("hTrackMatchedDEtaMCConversion",
1930        "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1931        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1932       fhTrackMatchedDEtaMCConversion[1]->SetYTitle("d#eta");
1933       fhTrackMatchedDEtaMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
1934       
1935       fhTrackMatchedDPhiMCConversion[1]  = new TH2F
1936       ("hTrackMatchedDPhiMCConversion",
1937        "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1938        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1939       fhTrackMatchedDPhiMCConversion[1]->SetYTitle("d#phi (rad)");
1940       fhTrackMatchedDPhiMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
1941       
1942       outputContainer->Add(fhTrackMatchedDEtaMCConversion[1]) ; 
1943       outputContainer->Add(fhTrackMatchedDPhiMCConversion[1]) ;
1944       
1945       
1946       fhTrackMatchedMCParticle[0]  = new TH2F
1947       ("hTrackMatchedMCParticleNoCut",
1948        "Origin of particle vs energy",
1949        nptbins,ptmin,ptmax,8,0,8); 
1950       fhTrackMatchedMCParticle[0]->SetXTitle("E (GeV)");   
1951       //fhTrackMatchedMCParticle[0]->SetYTitle("Particle type");
1952       
1953       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(1 ,"Photon");
1954       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(2 ,"Electron");
1955       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1956       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(4 ,"Rest");
1957       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1958       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1959       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1960       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1961       
1962       fhTrackMatchedMCParticle[1]  = new TH2F
1963       ("hTrackMatchedMCParticle",
1964        "Origin of particle vs energy",
1965        nptbins,ptmin,ptmax,8,0,8); 
1966       fhTrackMatchedMCParticle[1]->SetXTitle("E (GeV)");   
1967       //fhTrackMatchedMCParticle[1]->SetYTitle("Particle type");
1968       
1969       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(1 ,"Photon");
1970       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(2 ,"Electron");
1971       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1972       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(4 ,"Rest");
1973       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1974       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1975       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1976       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");      
1977       
1978       outputContainer->Add(fhTrackMatchedMCParticle[0]);            
1979       outputContainer->Add(fhTrackMatchedMCParticle[1]);       
1980       
1981     }
1982   }  
1983   
1984   if(fFillPileUpHistograms)
1985   {
1986     
1987     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1988     
1989     for(Int_t i = 0 ; i < 7 ; i++)
1990     {
1991       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
1992                                        Form("Cluster  p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1993       fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
1994       outputContainer->Add(fhPtPileUp[i]);
1995       
1996       fhPtChargedPileUp[i]  = new TH1F(Form("hPtChargedPileUp%s",pileUpName[i].Data()),
1997                                       Form("Charged clusters p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1998       fhPtChargedPileUp[i]->SetXTitle("p_{T} (GeV/c)");
1999       outputContainer->Add(fhPtChargedPileUp[i]);
2000
2001       fhPtPhotonPileUp[i]  = new TH1F(Form("hPtPhotonPileUp%s",pileUpName[i].Data()),
2002                                       Form("Selected photon p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2003       fhPtPhotonPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2004       outputContainer->Add(fhPtPhotonPileUp[i]);
2005       
2006       
2007       fhClusterEFracLongTimePileUp[i]  = new TH2F(Form("hClusterEFracLongTimePileUp%s",pileUpName[i].Data()),
2008                                              Form("Cluster E vs fraction of cluster energy from large T cells, %s Pile-Up event",pileUpName[i].Data()),
2009                                              nptbins,ptmin,ptmax,200,0,1);
2010       fhClusterEFracLongTimePileUp[i]->SetXTitle("E (GeV)");
2011       fhClusterEFracLongTimePileUp[i]->SetYTitle("E(large time) / E");
2012       outputContainer->Add(fhClusterEFracLongTimePileUp[i]);
2013       
2014       fhClusterTimeDiffPileUp[i]  = new TH2F(Form("hClusterTimeDiffPileUp%s",pileUpName[i].Data()),
2015                                 Form("Cluster E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2016                                              nptbins,ptmin,ptmax,200,-100,100);
2017       fhClusterTimeDiffPileUp[i]->SetXTitle("E (GeV)");
2018       fhClusterTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2019       outputContainer->Add(fhClusterTimeDiffPileUp[i]);
2020       
2021       fhClusterTimeDiffChargedPileUp[i]  = new TH2F(Form("hClusterTimeDiffChargedPileUp%s",pileUpName[i].Data()),
2022                                        Form("Charged clusters E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2023                                                     nptbins,ptmin,ptmax,200,-100,100);
2024       fhClusterTimeDiffChargedPileUp[i]->SetXTitle("E (GeV)");
2025       fhClusterTimeDiffChargedPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2026       outputContainer->Add(fhClusterTimeDiffChargedPileUp[i]);
2027       
2028       fhClusterTimeDiffPhotonPileUp[i]  = new TH2F(Form("hClusterTimeDiffPhotonPileUp%s",pileUpName[i].Data()),
2029                                       Form("Selected photon E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2030                                                    nptbins,ptmin,ptmax,200,-100,100);
2031       fhClusterTimeDiffPhotonPileUp[i]->SetXTitle("E (GeV)");
2032       fhClusterTimeDiffPhotonPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2033       outputContainer->Add(fhClusterTimeDiffPhotonPileUp[i]);
2034
2035       fhLambda0PileUp[i]  = new TH2F(Form("hLambda0PileUp%s",pileUpName[i].Data()),
2036                                      Form("Cluster E vs #lambda^{2}_{0} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2037                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2038       fhLambda0PileUp[i]->SetXTitle("E (GeV)");
2039       fhLambda0PileUp[i]->SetYTitle("#lambda^{2}_{0}");
2040       outputContainer->Add(fhLambda0PileUp[i]);
2041       
2042       fhLambda0ChargedPileUp[i]  = new TH2F(Form("hLambda0ChargedPileUp%s",pileUpName[i].Data()),
2043                                                     Form("Charged clusters E vs #lambda^{2}_{0}in cluster, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2044       fhLambda0ChargedPileUp[i]->SetXTitle("E (GeV)");
2045       fhLambda0ChargedPileUp[i]->SetYTitle("#lambda^{2}_{0}");
2046       outputContainer->Add(fhLambda0ChargedPileUp[i]);
2047
2048     }
2049     
2050     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2051     fhTimeENoCut->SetXTitle("E (GeV)");
2052     fhTimeENoCut->SetYTitle("time (ns)");
2053     outputContainer->Add(fhTimeENoCut);  
2054     
2055     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
2056     fhTimeESPD->SetXTitle("E (GeV)");
2057     fhTimeESPD->SetYTitle("time (ns)");
2058     outputContainer->Add(fhTimeESPD);  
2059     
2060     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
2061     fhTimeESPDMulti->SetXTitle("E (GeV)");
2062     fhTimeESPDMulti->SetYTitle("time (ns)");
2063     outputContainer->Add(fhTimeESPDMulti);  
2064     
2065     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50); 
2066     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2067     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
2068     outputContainer->Add(fhTimeNPileUpVertSPD);
2069     
2070     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 ); 
2071     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2072     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
2073     outputContainer->Add(fhTimeNPileUpVertTrack);  
2074     
2075     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50); 
2076     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2077     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
2078     outputContainer->Add(fhTimeNPileUpVertContributors);  
2079     
2080     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); 
2081     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
2082     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
2083     outputContainer->Add(fhTimePileUpMainVertexZDistance);  
2084     
2085     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50); 
2086     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
2087     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
2088     outputContainer->Add(fhTimePileUpMainVertexZDiamond); 
2089     
2090     TString title[] = {"no |t diff| cut","|t diff|<20 ns","|t diff|>20 ns","|t diff|>40 ns"};
2091     TString name [] = {"TDiffNoCut","TDiffSmaller20ns","TDiffLarger20ns","TDiffLarger40ns"};
2092     for(Int_t i = 0; i < 4; i++) 
2093     {      
2094       fhClusterMultSPDPileUp[i] = new TH2F(Form("fhClusterMultSPDPileUp_%s", name[i].Data()),
2095                                            Form("Number of clusters per pile up event with E > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
2096                                            nptbins,ptmin,ptmax,100,0,100); 
2097       fhClusterMultSPDPileUp[i]->SetYTitle("n clusters ");
2098       fhClusterMultSPDPileUp[i]->SetXTitle("E_{cluster max} (GeV)");
2099       outputContainer->Add(fhClusterMultSPDPileUp[i]) ;   
2100       
2101       fhClusterMultNoPileUp[i] = new TH2F(Form("fhClusterMultNoPileUp_%s", name[i].Data()),
2102                                           Form("Number of clusters per non pile up event with E > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
2103                                           nptbins,ptmin,ptmax,100,0,100); 
2104       fhClusterMultNoPileUp[i]->SetYTitle("n clusters ");
2105       fhClusterMultNoPileUp[i]->SetXTitle("E_{cluster max} (GeV)");
2106       outputContainer->Add(fhClusterMultNoPileUp[i]) ;         
2107     }
2108     
2109   }
2110   
2111   if(IsDataMC())
2112   {
2113     TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
2114                         "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
2115                         "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String"                                } ; 
2116     
2117     TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
2118                         "Conversion", "Hadron", "AntiNeutron","AntiProton",
2119                         "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
2120     
2121     for(Int_t i = 0; i < fNOriginHistograms; i++)
2122     { 
2123       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
2124                                 Form("cluster from %s : E ",ptype[i].Data()),
2125                                 nptbins,ptmin,ptmax); 
2126       fhMCE[i]->SetXTitle("E (GeV)");
2127       outputContainer->Add(fhMCE[i]) ; 
2128       
2129       fhMCPt[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
2130                            Form("cluster from %s : p_{T} ",ptype[i].Data()),
2131                            nptbins,ptmin,ptmax); 
2132       fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
2133       outputContainer->Add(fhMCPt[i]) ;
2134       
2135       fhMCEta[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
2136                            Form("cluster from %s : #eta ",ptype[i].Data()),
2137                            nptbins,ptmin,ptmax,netabins,etamin,etamax); 
2138       fhMCEta[i]->SetYTitle("#eta");
2139       fhMCEta[i]->SetXTitle("E (GeV)");
2140       outputContainer->Add(fhMCEta[i]) ;
2141       
2142       fhMCPhi[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
2143                            Form("cluster from %s : #phi ",ptype[i].Data()),
2144                            nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
2145       fhMCPhi[i]->SetYTitle("#phi (rad)");
2146       fhMCPhi[i]->SetXTitle("E (GeV)");
2147       outputContainer->Add(fhMCPhi[i]) ;
2148       
2149       
2150       fhMCDeltaE[i]  = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
2151                                  Form("MC - Reco E from %s",pname[i].Data()), 
2152                                  nptbins,ptmin,ptmax, 200,-50,50); 
2153       fhMCDeltaE[i]->SetXTitle("#Delta E (GeV)");
2154       outputContainer->Add(fhMCDeltaE[i]);
2155       
2156       fhMCDeltaPt[i]  = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
2157                                   Form("MC - Reco p_{T} from %s",pname[i].Data()), 
2158                                   nptbins,ptmin,ptmax, 200,-50,50); 
2159       fhMCDeltaPt[i]->SetXTitle("#Delta p_{T} (GeV/c)");
2160       outputContainer->Add(fhMCDeltaPt[i]);
2161             
2162       fhMC2E[i]  = new TH2F (Form("h2E_MC%s",pname[i].Data()),
2163                              Form("E distribution, reconstructed vs generated from %s",pname[i].Data()), 
2164                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
2165       fhMC2E[i]->SetXTitle("E_{rec} (GeV)");
2166       fhMC2E[i]->SetYTitle("E_{gen} (GeV)");
2167       outputContainer->Add(fhMC2E[i]);          
2168       
2169       fhMC2Pt[i]  = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
2170                               Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()), 
2171                               nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
2172       fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/c)");
2173       fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/c)");
2174       outputContainer->Add(fhMC2Pt[i]);
2175       
2176       
2177     }
2178     
2179     TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
2180                          "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ; 
2181     
2182     TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
2183                          "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
2184     
2185     for(Int_t i = 0; i < fNPrimaryHistograms; i++)
2186     { 
2187       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
2188                            Form("primary photon %s : E ",pptype[i].Data()),
2189                            nptbins,ptmin,ptmax); 
2190       fhEPrimMC[i]->SetXTitle("E (GeV)");
2191       outputContainer->Add(fhEPrimMC[i]) ; 
2192       
2193       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
2194                             Form("primary photon %s : p_{T} ",pptype[i].Data()),
2195                             nptbins,ptmin,ptmax); 
2196       fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
2197       outputContainer->Add(fhPtPrimMC[i]) ;
2198       
2199       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
2200                              Form("primary photon %s : Rapidity ",pptype[i].Data()),
2201                              nptbins,ptmin,ptmax,800,-8,8); 
2202       fhYPrimMC[i]->SetYTitle("Rapidity");
2203       fhYPrimMC[i]->SetXTitle("E (GeV)");
2204       outputContainer->Add(fhYPrimMC[i]) ;
2205       
2206       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
2207                              Form("primary photon %s : #phi ",pptype[i].Data()),
2208                              nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
2209       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
2210       fhPhiPrimMC[i]->SetXTitle("E (GeV)");
2211       outputContainer->Add(fhPhiPrimMC[i]) ;
2212      
2213       
2214       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
2215                                Form("primary photon %s in acceptance: E ",pptype[i].Data()),
2216                                nptbins,ptmin,ptmax); 
2217       fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
2218       outputContainer->Add(fhEPrimMCAcc[i]) ; 
2219       
2220       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
2221                                 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
2222                                 nptbins,ptmin,ptmax); 
2223       fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
2224       outputContainer->Add(fhPtPrimMCAcc[i]) ;
2225       
2226       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
2227                                  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
2228                                  nptbins,ptmin,ptmax,100,-1,1); 
2229       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
2230       fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
2231       outputContainer->Add(fhYPrimMCAcc[i]) ;
2232       
2233       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
2234                                  Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
2235                                  nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
2236       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
2237       fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
2238       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
2239       
2240     }
2241                 
2242     if(fFillSSHistograms)
2243     {
2244       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ; 
2245       
2246       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
2247       
2248       for(Int_t i = 0; i < 6; i++)
2249       { 
2250         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
2251                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
2252                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2253         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
2254         fhMCELambda0[i]->SetXTitle("E (GeV)");
2255         outputContainer->Add(fhMCELambda0[i]) ; 
2256         
2257         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
2258                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
2259                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2260         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
2261         fhMCELambda1[i]->SetXTitle("E (GeV)");
2262         outputContainer->Add(fhMCELambda1[i]) ; 
2263         
2264         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
2265                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
2266                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2267         fhMCEDispersion[i]->SetYTitle("D^{2}");
2268         fhMCEDispersion[i]->SetXTitle("E (GeV)");
2269         outputContainer->Add(fhMCEDispersion[i]) ; 
2270         
2271         fhMCNCellsE[i]  = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
2272                                     Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()), 
2273                                     nptbins,ptmin,ptmax, nbins,nmin,nmax); 
2274         fhMCNCellsE[i]->SetXTitle("E (GeV)");
2275         fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
2276         outputContainer->Add(fhMCNCellsE[i]);  
2277         
2278         fhMCMaxCellDiffClusterE[i]  = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
2279                                                 Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
2280                                                 nptbins,ptmin,ptmax, 500,0,1.); 
2281         fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
2282         fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2283         outputContainer->Add(fhMCMaxCellDiffClusterE[i]);  
2284         
2285         if(!fFillOnlySimpleSSHisto)
2286         {
2287           fhMCLambda0vsClusterMaxCellDiffE0[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
2288                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
2289                                                            ssbins,ssmin,ssmax,500,0,1.); 
2290           fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
2291           fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2292           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ; 
2293           
2294           fhMCLambda0vsClusterMaxCellDiffE2[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
2295                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
2296                                                            ssbins,ssmin,ssmax,500,0,1.); 
2297           fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
2298           fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2299           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ; 
2300           
2301           fhMCLambda0vsClusterMaxCellDiffE6[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
2302                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
2303                                                            ssbins,ssmin,ssmax,500,0,1.); 
2304           fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
2305           fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2306           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ; 
2307           
2308           fhMCNCellsvsClusterMaxCellDiffE0[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
2309                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
2310                                                           nbins/5,nmin,nmax/5,500,0,1.); 
2311           fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
2312           fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2313           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ; 
2314           
2315           fhMCNCellsvsClusterMaxCellDiffE2[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
2316                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
2317                                                           nbins/5,nmin,nmax/5,500,0,1.); 
2318           fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
2319           fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2320           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ; 
2321           
2322           fhMCNCellsvsClusterMaxCellDiffE6[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
2323                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
2324                                                           nbins/5,nmin,nmax/5,500,0,1.); 
2325           fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
2326           fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
2327           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ; 
2328           
2329           if(fCalorimeter=="EMCAL")
2330           {
2331             fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
2332                                          Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
2333                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
2334             fhMCEDispEta[i]->SetXTitle("E (GeV)");
2335             fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
2336             outputContainer->Add(fhMCEDispEta[i]);     
2337             
2338             fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
2339                                          Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
2340                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
2341             fhMCEDispPhi[i]->SetXTitle("E (GeV)");
2342             fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2343             outputContainer->Add(fhMCEDispPhi[i]);  
2344             
2345             fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
2346                                            Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data()),  
2347                                            nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
2348             fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
2349             fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
2350             outputContainer->Add(fhMCESumEtaPhi[i]);
2351             
2352             fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
2353                                                 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),  
2354                                                 nptbins,ptmin,ptmax,200,-10,10); 
2355             fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
2356             fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
2357             outputContainer->Add(fhMCEDispEtaPhiDiff[i]);    
2358             
2359             fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
2360                                             Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data()),  
2361                                             nptbins,ptmin,ptmax, 200,-1,1); 
2362             fhMCESphericity[i]->SetXTitle("E (GeV)");
2363             fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
2364             outputContainer->Add(fhMCESphericity[i]);
2365             
2366             for(Int_t ie = 0; ie < 7; ie++)
2367             {
2368               fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
2369                                                     Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]), 
2370                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
2371               fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2372               fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2373               outputContainer->Add(fhMCDispEtaDispPhi[ie][i]); 
2374               
2375               fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pnamess[i].Data()),
2376                                                     Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]), 
2377                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
2378               fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
2379               fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2380               outputContainer->Add(fhMCLambda0DispEta[ie][i]);       
2381               
2382               fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
2383                                                     Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]), 
2384                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
2385               fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
2386               fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2387               outputContainer->Add(fhMCLambda0DispPhi[ie][i]); 
2388             }
2389           }
2390         }
2391       }// loop    
2392       
2393       if(!GetReader()->IsEmbeddedClusterSelectionOn())
2394       {
2395         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
2396                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
2397                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2398         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
2399         fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
2400         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ; 
2401         
2402         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
2403                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
2404                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2405         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
2406         fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
2407         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ; 
2408         
2409         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
2410                                                "cluster from Photon : E vs #lambda_{0}^{2}",
2411                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2412         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
2413         fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
2414         outputContainer->Add(fhMCPhotonELambda0NOverlap) ; 
2415         
2416       } //No embedding     
2417       
2418       //Fill histograms to check shape of embedded clusters
2419       if(GetReader()->IsEmbeddedClusterSelectionOn())
2420       {
2421         
2422         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
2423                                                    "Energy Fraction of embedded signal versus cluster energy",
2424                                                    nptbins,ptmin,ptmax,100,0.,1.); 
2425         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
2426         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
2427         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ; 
2428         
2429         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
2430                                                     "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
2431                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2432         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
2433         fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
2434         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ; 
2435         
2436         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
2437                                                       "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
2438                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2439         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
2440         fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
2441         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ; 
2442         
2443         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
2444                                                    "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
2445                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2446         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2447         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
2448         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ; 
2449         
2450         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
2451                                                  "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2452                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2453         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2454         fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
2455         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ; 
2456         
2457         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
2458                                                  "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
2459                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2460         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
2461         fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
2462         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ; 
2463         
2464         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
2465                                                    "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
2466                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2467         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
2468         fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
2469         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ; 
2470         
2471         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
2472                                                 "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
2473                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2474         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2475         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
2476         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ; 
2477         
2478         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
2479                                               "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2480                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2481         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2482         fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
2483         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ; 
2484         
2485       }// embedded histograms
2486       
2487       
2488     }// Fill SS MC histograms
2489     
2490   }//Histos with MC
2491       
2492   return outputContainer ;
2493   
2494 }
2495
2496 //_______________________
2497 void AliAnaPhoton::Init()
2498 {
2499   
2500   //Init
2501   //Do some checks
2502   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
2503   {
2504     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2505     abort();
2506   }
2507   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
2508   {
2509     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2510     abort();
2511   }
2512   
2513   if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
2514   
2515 }
2516
2517 //____________________________________________________________________________
2518 void AliAnaPhoton::InitParameters()
2519 {
2520   
2521   //Initialize the parameters of the analysis.
2522   AddToHistogramsName("AnaPhoton_");
2523   
2524   fCalorimeter = "EMCAL" ;
2525   fMinDist     = 2.;
2526   fMinDist2    = 4.;
2527   fMinDist3    = 5.;
2528         
2529   fTimeCutMin  =-1000000;
2530   fTimeCutMax  = 1000000;
2531   fNCellsCut   = 0;
2532         
2533   fRejectTrackMatch       = kTRUE ;
2534         
2535 }
2536
2537 //__________________________________________________________________
2538 void  AliAnaPhoton::MakeAnalysisFillAOD() 
2539 {
2540   //Do photon analysis and fill aods
2541   
2542   //Get the vertex 
2543   Double_t v[3] = {0,0,0}; //vertex ;
2544   GetReader()->GetVertex(v);
2545   
2546   //Select the Calorimeter of the photon
2547   TObjArray * pl = 0x0; 
2548   AliVCaloCells* cells    = 0;  
2549   if      (fCalorimeter == "PHOS" )
2550   {
2551     pl    = GetPHOSClusters();
2552     cells = GetPHOSCells();
2553   }
2554   else if (fCalorimeter == "EMCAL")
2555   {
2556     pl    = GetEMCALClusters();
2557     cells = GetEMCALCells();
2558   }
2559   
2560   if(!pl) 
2561   {
2562     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2563     return;
2564   }
2565   
2566   FillPileUpHistogramsPerEvent(pl); 
2567   
2568   // Loop on raw clusters before filtering in the reader and fill control histogram
2569   if((GetReader()->GetEMCALClusterListName()=="" && fCalorimeter=="EMCAL") || fCalorimeter=="PHOS")
2570   {
2571     for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
2572     {
2573       AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
2574       if     (fCalorimeter == "PHOS"  && clus->IsPHOS()  && clus->E() > GetReader()->GetPHOSPtMin() ) fhClusterCuts[0]->Fill(clus->E());
2575       else if(fCalorimeter == "EMCAL" && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin()) fhClusterCuts[0]->Fill(clus->E());
2576     }
2577   }
2578   else 
2579   { // reclusterized
2580     TClonesArray * clusterList = 0;
2581
2582     if(GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()))
2583       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
2584     else if(GetReader()->GetOutputEvent())
2585       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
2586     
2587     if(clusterList)
2588     {
2589       Int_t nclusters = clusterList->GetEntriesFast();
2590       for (Int_t iclus =  0; iclus <  nclusters; iclus++) 
2591       {
2592         AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));        
2593         if(clus)fhClusterCuts[0]->Fill(clus->E());
2594       }  
2595     }
2596   }
2597   
2598   //Init arrays, variables, get number of clusters
2599   TLorentzVector mom, mom2 ;
2600   Int_t nCaloClusters = pl->GetEntriesFast();
2601   
2602   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
2603   
2604   //----------------------------------------------------
2605   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2606   //----------------------------------------------------
2607   // Loop on clusters
2608   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
2609   {    
2610           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo)); 
2611     //printf("calo %d, %f\n",icalo,calo->E());
2612     
2613     //Get the index where the cluster comes, to retrieve the corresponding vertex
2614     Int_t evtIndex = 0 ; 
2615     if (GetMixedEvent()) 
2616     {
2617       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
2618       //Get the vertex and check it is not too large in z
2619       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
2620     }
2621     
2622     //Cluster selection, not charged, with photon id and in fiducial cut          
2623     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2624     {
2625       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
2626     else
2627     {
2628       Double_t vertex[]={0,0,0};
2629       calo->GetMomentum(mom,vertex) ;
2630     }
2631         
2632     //--------------------------------------
2633     // Cluster selection
2634     //--------------------------------------
2635     Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
2636     if(!ClusterSelected(calo,mom,nMaxima)) continue;
2637
2638     //----------------------------
2639     //Create AOD for analysis
2640     //----------------------------
2641     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
2642     
2643     //...............................................
2644     //Set the indeces of the original caloclusters (MC, ID), and calorimeter  
2645     Int_t label = calo->GetLabel();
2646     aodph.SetLabel(label);
2647     aodph.SetCaloLabel(calo->GetID(),-1);
2648     aodph.SetDetector(fCalorimeter);
2649     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
2650     
2651     //...............................................
2652     //Set bad channel distance bit
2653     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2654     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
2655     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; 
2656     else                         aodph.SetDistToBad(0) ;
2657     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
2658     
2659     //--------------------------------------------------------------------------------------
2660     // Play with the MC stack if available
2661     //--------------------------------------------------------------------------------------
2662     
2663     //Check origin of the candidates
2664     Int_t tag = -1;
2665     
2666     if(IsDataMC())
2667     {
2668       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex());
2669       aodph.SetTag(tag);
2670       
2671       if(GetDebug() > 0)
2672         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
2673     }//Work with stack also   
2674         
2675     
2676     //--------------------------------------------------------------------------------------
2677     //Fill some shower shape histograms before PID is applied
2678     //--------------------------------------------------------------------------------------
2679     
2680     FillShowerShapeHistograms(calo,tag);
2681     
2682     //-------------------------------------
2683     //PID selection or bit setting
2684     //-------------------------------------
2685     
2686     //...............................................
2687     // Data, PID check on
2688     if(IsCaloPIDOn())
2689     {
2690       // Get most probable PID, 2 options check bayesian PID weights or redo PID
2691       // By default, redo PID
2692       
2693       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));
2694       
2695       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
2696       
2697       //If cluster does not pass pid, not photon, skip it.
2698       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;                   
2699       
2700     }
2701
2702     //...............................................
2703     // Data, PID check off
2704     else
2705     {
2706       //Set PID bits for later selection (AliAnaPi0 for example)
2707       //GetIdentifiedParticleType already called in SetPIDBits.
2708       
2709       GetCaloPID()->SetPIDBits(calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
2710       
2711       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");               
2712     }
2713     
2714     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
2715                               aodph.Pt(), aodph.GetIdentifiedParticleType());
2716     
2717     fhClusterCuts[9]->Fill(calo->E());
2718
2719     fhNLocMax->Fill(calo->E(),nMaxima);
2720
2721     // Matching after cuts
2722     if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,1);
2723     
2724     // Fill histograms to undertand pile-up before other cuts applied
2725     // Remember to relax time cuts in the reader
2726     FillPileUpHistograms(calo->E(),mom.Pt(),calo->GetTOF()*1e9);
2727     
2728     // Add number of local maxima to AOD, method name in AOD to be FIXED
2729     aodph.SetFiducialArea(nMaxima);
2730     
2731
2732     //Add AOD with photon object to aod branch
2733     AddAODParticle(aodph);
2734     
2735   }//loop
2736         
2737   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());  
2738   
2739 }
2740
2741 //__________________________________________________________________
2742 void  AliAnaPhoton::MakeAnalysisFillHistograms() 
2743 {
2744   //Fill histograms
2745       
2746   // Get vertex
2747   Double_t v[3] = {0,0,0}; //vertex ;
2748   GetReader()->GetVertex(v);
2749   //fhVertex->Fill(v[0],v[1],v[2]);  
2750   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2751   
2752   //----------------------------------
2753   //Loop on stored AOD photons
2754   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2755   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2756   
2757   for(Int_t iaod = 0; iaod < naod ; iaod++)
2758   {
2759     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2760     Int_t pdg = ph->GetIdentifiedParticleType();
2761     
2762     if(GetDebug() > 3) 
2763       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", 
2764              ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
2765     
2766     //If PID used, fill histos with photons in Calorimeter fCalorimeter
2767     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue; 
2768     if(ph->GetDetector() != fCalorimeter) continue;
2769     
2770     if(GetDebug() > 2) 
2771       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
2772     
2773     //................................
2774     //Fill photon histograms 
2775     Float_t ptcluster  = ph->Pt();
2776     Float_t phicluster = ph->Phi();
2777     Float_t etacluster = ph->Eta();
2778     Float_t ecluster   = ph->E();
2779     
2780     fhEPhoton   ->Fill(ecluster);
2781     fhPtPhoton  ->Fill(ptcluster);
2782     fhPhiPhoton ->Fill(ptcluster,phicluster);
2783     fhEtaPhoton ->Fill(ptcluster,etacluster);    
2784     if     (ecluster   > 0.5) fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
2785     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
2786
2787     //Get original cluster, to recover some information
2788     Int_t absID             = 0; 
2789     Float_t maxCellFraction = 0;
2790     AliVCaloCells* cells    = 0;
2791     TObjArray * clusters    = 0; 
2792     if(fCalorimeter == "EMCAL")
2793     {
2794       cells    = GetEMCALCells();
2795       clusters = GetEMCALClusters();
2796     }
2797     else
2798     {
2799       cells    = GetPHOSCells();
2800       clusters = GetPHOSClusters();
2801     }
2802     
2803     Int_t iclus = -1;
2804     AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus); 
2805     if(cluster)
2806     {
2807       absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
2808       
2809       // Control histograms
2810       fhMaxCellDiffClusterE->Fill(ph->E(),maxCellFraction);
2811       fhNCellsE            ->Fill(ph->E(),cluster->GetNCells());
2812       fhTimeE              ->Fill(ph->E(),cluster->GetTOF()*1.e9);
2813       if(cells)
2814       {
2815       for(Int_t icell = 0; icell <  cluster->GetNCells(); icell++)
2816         fhCellsE->Fill(ph->E(),cells->GetCellAmplitude(cluster->GetCellsAbsId()[icell]));
2817       }
2818     }
2819     
2820     //.......................................
2821     //Play with the MC data if available
2822     if(IsDataMC())
2823     {
2824       if(GetDebug()>0)
2825       {
2826         if(GetReader()->ReadStack() && !GetMCStack())
2827         {
2828           printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called?\n");
2829         }
2830         else if(GetReader()->ReadAODMCParticles() && !GetReader()->GetAODMCParticles(0))
2831         {
2832           printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
2833         }       
2834       }    
2835       
2836       FillAcceptanceHistograms();
2837       
2838       //....................................................................
2839       // Access MC information in stack if requested, check that it exists.
2840       Int_t label =ph->GetLabel();
2841       
2842       if(label < 0) 
2843       {
2844         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
2845         continue;
2846       }
2847       
2848       Float_t eprim   = 0;
2849       Float_t ptprim  = 0;
2850       Bool_t ok = kFALSE;
2851       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2852       if(ok)
2853       {
2854         eprim   = primary.Energy();
2855         ptprim  = primary.Pt();         
2856       }
2857       
2858       Int_t tag =ph->GetTag();
2859       Int_t mcParticleTag = -1;
2860       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[kmcPhoton])
2861       {
2862         fhMCE  [kmcPhoton] ->Fill(ecluster);
2863         fhMCPt [kmcPhoton] ->Fill(ptcluster);
2864         fhMCPhi[kmcPhoton] ->Fill(ecluster,phicluster);
2865         fhMCEta[kmcPhoton] ->Fill(ecluster,etacluster);
2866         
2867         fhMC2E     [kmcPhoton] ->Fill(ecluster, eprim);
2868         fhMC2Pt    [kmcPhoton] ->Fill(ptcluster, ptprim);     
2869         fhMCDeltaE [kmcPhoton] ->Fill(ecluster,eprim-ecluster);
2870         fhMCDeltaPt[kmcPhoton] ->Fill(ptcluster,ptprim-ptcluster); 
2871         
2872         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && 
2873            GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)     &&
2874            fhMCE[kmcConversion])
2875         {
2876           fhMCE  [kmcConversion] ->Fill(ecluster);
2877           fhMCPt [kmcConversion] ->Fill(ptcluster);
2878           fhMCPhi[kmcConversion] ->Fill(ecluster,phicluster);
2879           fhMCEta[kmcConversion] ->Fill(ecluster,etacluster);
2880           
2881           fhMC2E     [kmcConversion] ->Fill(ecluster, eprim);
2882           fhMC2Pt    [kmcConversion] ->Fill(ptcluster, ptprim);     
2883           fhMCDeltaE [kmcConversion] ->Fill(ecluster,eprim-ecluster);
2884           fhMCDeltaPt[kmcConversion] ->Fill(ptcluster,ptprim-ptcluster); 
2885         }                       
2886         
2887         if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[kmcPrompt])
2888         {
2889           mcParticleTag = kmcPrompt;
2890         }
2891         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[kmcFragmentation])
2892         {
2893           mcParticleTag = kmcFragmentation;
2894         }
2895         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[kmcISR])
2896         {
2897           mcParticleTag = kmcISR; 
2898         }
2899         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) && 
2900                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[kmcPi0Decay])
2901         {
2902           mcParticleTag = kmcPi0Decay;
2903         }
2904         else if((( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) && 
2905                   !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)        ) || 
2906                    GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[kmcOtherDecay])
2907         {
2908           mcParticleTag = kmcOtherDecay;
2909         }
2910         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[kmcPi0])
2911         {
2912           mcParticleTag = kmcPi0;   
2913         } 
2914         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[kmcEta])
2915         {
2916           mcParticleTag = kmcEta;
2917         }      
2918       }
2919       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[kmcAntiNeutron])
2920       {
2921         mcParticleTag = kmcAntiNeutron;
2922       }
2923       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[kmcAntiProton])
2924       {
2925         mcParticleTag = kmcAntiProton; 
2926       } 
2927       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[kmcElectron])
2928       {
2929         mcParticleTag = kmcElectron;            
2930       }     
2931       else if( fhMCE[kmcOther])
2932       {
2933         mcParticleTag = kmcOther;
2934         
2935         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2936         //                                      ph->GetLabel(),ph->Pt());
2937         //                for(Int_t i = 0; i < 20; i++) {
2938         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2939         //                }
2940         //                printf("\n");
2941         
2942       }
2943       
2944       fhMCE  [mcParticleTag] ->Fill(ecluster);
2945       fhMCPt [mcParticleTag] ->Fill(ptcluster);
2946       fhMCPhi[mcParticleTag] ->Fill(ecluster,phicluster);
2947       fhMCEta[mcParticleTag] ->Fill(ecluster,etacluster);
2948       
2949       fhMC2E[mcParticleTag]     ->Fill(ecluster, eprim);
2950       fhMC2Pt[mcParticleTag]    ->Fill(ptcluster, ptprim);     
2951       fhMCDeltaE[mcParticleTag] ->Fill(ecluster,eprim-ecluster);
2952       fhMCDeltaPt[mcParticleTag]->Fill(ptcluster,ptprim-ptcluster); 
2953       
2954     }//Histograms with MC
2955     
2956   }// aod loop
2957   
2958 }
2959
2960
2961 //__________________________________________________________________
2962 void AliAnaPhoton::Print(const Option_t * opt) const
2963 {
2964   //Print some relevant parameters set for the analysis
2965   
2966   if(! opt)
2967     return;
2968   
2969   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2970   AliAnaCaloTrackCorrBaseClass::Print(" ");
2971
2972   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
2973   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
2974   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2975   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2976   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
2977   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
2978   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
2979   printf("    \n") ;
2980         
2981