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