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