AliPhysics  608b256 (608b256)
AliAnaPhoton.cxx
Go to the documentation of this file.
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 // --- ROOT system ---
17 #include <TH3F.h>
18 #include <TH2F.h>
19 #include <TClonesArray.h>
20 #include <TObjString.h>
21 #include "AliVParticle.h"
22 #include "TDatabasePDG.h"
23 
24 // --- Analysis system ---
25 #include "AliAnaPhoton.h"
26 #include "AliCaloTrackReader.h"
27 #include "AliMCEvent.h"
28 #include "AliCaloPID.h"
29 #include "AliMCAnalysisUtils.h"
30 #include "AliFiducialCut.h"
31 #include "AliVCluster.h"
32 #include "AliMixedEvent.h"
33 #include "AliAODEvent.h"
34 #include "AliESDEvent.h"
35 #include "AliMCEvent.h"
36 
37 // --- Detectors ---
38 #include "AliPHOSGeoUtils.h"
39 #include "AliEMCALGeometry.h"
40 
42 ClassImp(AliAnaPhoton) ;
44 
45 //____________________________
48 //____________________________
51 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
52 fRejectTrackMatch(0), fFillTMHisto(kFALSE), fFillTMHistoTrackPt(0),
53 fTimeCutMin(-10000), fTimeCutMax(10000),
54 fNCellsCut(0),
55 fNLMCutMin(-1), fNLMCutMax(10),
56 fFillSSHistograms(0), fFillSSPerSMHistograms(0), fFillEMCALRegionSSHistograms(0),
57 fFillConversionVertexHisto(0),fFillOnlySimpleSSHisto(1),
58 fFillSSNLocMaxHisto(0),
59 fFillTrackMultHistograms(0),
60 fNOriginHistograms(9), fNPrimaryHistograms(5),
61 fMomentum(), fMomentum2(),
62 fPrimaryMom(), fProdVertex(),
63 fConstantTimeShift(0), fFillEBinAcceptanceHisto(0), fNEBinCuts(0),
64 fStudyActivityNearCluster(0),
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), fhDispPt(0),
79 fhLam0E(0), fhLam0Pt(0),
80 fhLam1E(0), fhLam1Pt(0),
81 fhLam0PtNLM1(0), fhLam0PtNLM2(0),
82 fhLam1PtNLM1(0), fhLam1PtNLM2(0),
83 fhDispETRD(0), fhLam0ETRD(0), fhLam0PtTRD(0), fhLam1ETRD(0),
84 fhDispETM(0), fhLam0ETM(0), fhLam0PtTM(0), fhLam1ETM(0),
85 fhDispETMTRD(0), fhLam0ETMTRD(0), fhLam0PtTMTRD(0), fhLam1ETMTRD(0),
86 
87 fhNCellsLam0LowE(0), fhNCellsLam1LowE(0), fhNCellsDispLowE(0),
88 fhNCellsLam0HighE(0), fhNCellsLam1HighE(0), fhNCellsDispHighE(0),
89 
90 fhEtaLam0LowE(0), fhPhiLam0LowE(0),
91 fhEtaLam0HighE(0), fhPhiLam0HighE(0),
92 fhLam0DispLowE(0), fhLam0DispHighE(0),
93 fhLam1Lam0LowE(0), fhLam1Lam0HighE(0),
94 fhDispLam1LowE(0), fhDispLam1HighE(0),
95 fhDispEtaE(0), fhDispPhiE(0),
96 fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
97 fhDispEtaPhiDiffE(0), fhSphericityE(0),
98 fhDispSumEtaDiffE(0), fhDispSumPhiDiffE(0),
99 
100 // MC histograms
101 //fhMCPhotonELambda0NoOverlap(0), fhMCPhotonELambda0TwoOverlap(0), fhMCPhotonELambda0NOverlap(0),
102 // Embedding
103 fhEmbeddedSignalFractionEnergy(0),
104 fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonELambda0MostlySignal(0),
105 fhEmbedPhotonELambda0MostlyBkg(0), fhEmbedPhotonELambda0FullBkg(0),
106 fhEmbedPi0ELambda0FullSignal(0), fhEmbedPi0ELambda0MostlySignal(0),
107 fhEmbedPi0ELambda0MostlyBkg(0), fhEmbedPi0ELambda0FullBkg(0),
108 
109 fhTimePtPhotonNoCut(0), fhTimePtPhotonSPD(0),
110 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
111 fhPtPhotonNPileUpSPDVtx(0), fhPtPhotonNPileUpTrkVtx(0),
112 fhPtPhotonNPileUpSPDVtxTimeCut(0), fhPtPhotonNPileUpTrkVtxTimeCut(0),
113 fhPtPhotonNPileUpSPDVtxTimeCut2(0), fhPtPhotonNPileUpTrkVtxTimeCut2(0),
114 
115 fhEClusterSM(0), fhEPhotonSM(0),
116 fhPtClusterSM(0), fhPtPhotonSM(0),
117 fhMCConversionVertex(0), fhMCConversionVertexTRD(0),
118 //fhDistanceAddedPhotonAddedPrimarySignal (0), fhDistanceHijingPhotonAddedPrimarySignal (0),
119 //fhDistanceAddedPhotonAddedSecondarySignal(0), fhDistanceHijingPhotonAddedSecondarySignal(0),
120 //fhDistanceAddedPhotonHijingSecondary(0)
121 fhLocalRegionClusterEnergySumHijing2(0),
122 fhLocalRegionClusterMultiplicityHijing2(0),
123 fhLocalRegionClusterEnergySumPerCentralityHijing2(0),
124 fhLocalRegionClusterMultiplicityPerCentralityHijing2(0),
125 fhDistance2AddedSignals(0), fhDistanceAddedSignalsHijing(0),
126 fhDistance2Hijing(0)
127 {
128  for(Int_t i = 0; i < fgkNmcTypes; i++)
129  {
130  fhMCPt [i] = 0;
131  fhMCE [i] = 0;
132  fhMCPhi [i] = 0;
133  fhMCEta [i] = 0;
134  fhMCDeltaE [i] = 0;
135  fhMCDeltaPt[i] = 0;
136  fhMC2E [i] = 0;
137  fhMC2Pt [i] = 0;
138  }
139 
140  for(Int_t i = 0; i < fgkNmcPrimTypes; i++)
141  {
142  fhPtPrimMC [i] = 0;
143  fhEPrimMC [i] = 0;
144  fhPhiPrimMC[i] = 0;
145  fhEtaPrimMC[i] = 0;
146  fhYPrimMC [i] = 0;
147 
148  fhPtPrimMCAcc [i] = 0;
149  fhEPrimMCAcc [i] = 0;
150  fhPhiPrimMCAcc[i] = 0;
151  fhEtaPrimMCAcc[i] = 0;
152  fhYPrimMCAcc [i] = 0;
153  }
154 
155  for(Int_t i = 0; i < 7; i++)
156  {
157  fhDispEtaDispPhi[i] = 0;
158  fhLambda0DispPhi[i] = 0;
159  fhLambda0DispEta[i] = 0;
160 
161  fhPtPhotonPileUp[i] = 0;
163 
164  for(Int_t j = 0; j < fgkNssTypes; j++)
165  {
166  fhMCDispEtaDispPhi[i][j] = 0;
167  fhMCLambda0DispEta[i][j] = 0;
168  fhMCLambda0DispPhi[i][j] = 0;
169  }
170  }
171 
172  for(Int_t i = 0; i < fgkNssTypes; i++)
173  {
174  fhMCELambda0 [i] = 0;
175  fhMCPtLambda0 [i] = 0;
176  fhMCELambda1 [i] = 0;
177  fhMCEDispersion [i] = 0;
178  fhMCNCellsE [i] = 0;
180  fhLambda0DispEta[i] = 0;
181  fhLambda0DispPhi[i] = 0;
182 
183  for(Int_t iover = 0 ; iover < 3; iover++)
184  fhMCPtLambda0Overlaps[i][iover] = 0;
185 
192 
193  fhMCEDispEta [i] = 0;
194  fhMCEDispPhi [i] = 0;
195  fhMCESumEtaPhi [i] = 0;
196  fhMCEDispEtaPhiDiff[i] = 0;
197  fhMCESphericity [i] = 0;
198  }
199 
200  for(Int_t i = 0; i < 5; i++)
201  {
202  fhClusterCutsE [i] = 0;
203  fhClusterCutsPt[i] = 0;
204  }
205 
206  for(Int_t icut = 0; icut < 10; icut++)
207  {
208  fhPtPhotonNTracks [icut] = 0;
209  fhPtPhotonSumPtTracks[icut] = 0;
210  }
211 
212  // Track matching residuals
213  for(Int_t i = 0; i < 2; i++)
214  {
215  //fhTrackMatchedDEta [i] = 0; fhTrackMatchedDPhi [i] = 0; fhTrackMatchedDEtaDPhi [i] = 0;
218 
219  //fhTrackMatchedDEtaTrackPt [i] = 0; fhTrackMatchedDPhiTrackPt [i] = 0; fhTrackMatchedDEtaDPhiTrackPt [i] = 0;
222 
232  fhdEdx[i] = 0; fhEOverP[i] = 0;
233  fhdEdxTrackPt[i] = 0; fhEOverPTrackPt[i] = 0;
234  fhEOverPTRD[i] = 0;
235  }
236 
237  for(Int_t i = 0; i < 4; i++)
238  {
239  fhMCParticle [i] = 0;
240  fhMCParticleConverted [i] = 0;
241  fhMCParticleVsErecEgen [i] = 0;
242  fhMCParticleVsNOverlaps[i] = 0;
243  }
244 
245  for(Int_t i = 0; i < 6; i++)
246  {
247  fhMCConversionLambda0Rcut [i] = 0;
249  fhMCConversionLambda1Rcut [i] = 0;
251  }
252 
253  for(Int_t ieta = 0; ieta < 4; ieta++)
254  {
255  for(Int_t iphi = 0; iphi < 3; iphi++)
256  {
257 // fhLam0EMCALRegion [ieta][iphi] = 0;
258 // fhLam0EMCALRegionTRD[ieta][iphi] = 0;
259 //
260 // for(Int_t i = 0; i < 6; i++)
261 // {
262 // fhLam0EMCALRegionMCConvRcut [ieta][iphi][i] = 0;
263 // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][i] = 0;
264 // }
265  for(Int_t ism =0; ism < 20; ism++)
266  {
267  fhLam0EMCALRegionPerSM[ieta][iphi][ism] = 0;
268  fhLam1EMCALRegionPerSM[ieta][iphi][ism] = 0;
269  }
270  }
271  }
272 
273  for(Int_t il0 = 0; il0 < 2; il0++)
274  {
275  for(Int_t i = 0; i < 7; i++)
276  {
277  fhEtaPhiLam0BinPtBin [il0][i] = 0 ;
278  fhEtaPhiLargeTimeInClusterCell[il0][i] = 0 ;
279 
280  fhColRowLam0BinPtBin [il0][i] = 0 ;
281  fhColRowLam0BinPtBinWeighted [il0][i] = 0 ;
282  fhColRowLam0BinPtBinLargeTime [il0][i] = 0 ;
283  fhCellClusterIndexEAndTime [il0][i] = 0 ;
284  fhCellClusterEAndTime [il0][i] = 0 ;
285  fhCellClusterEFracAndTime [il0][i] = 0 ;
286  }
287 
288  for(Int_t ism =0; ism < 20; ism++)
289  {
290  fhLam1Lam0BinPerSM [il0][ism] = 0;
291  fhTimeLam0BinPerSM [il0][ism] = 0;
292  fhTimeLam0BinPerSMWeighted [il0][ism] = 0;
293  fhDTimeLam0BinPerSM [il0][ism] = 0;
294  fhDTimeLam0BinPerSMWeighted [il0][ism] = 0;
295  fhCellClusterEFracLam0BinPerSM [il0][ism] = 0;
296 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism] = 0;
297  fhCellClusterELam0BinPerSM [il0][ism] = 0;
299 
302  fhCellClusterELam0BinPerSMLargeTime [il0][ism] = 0;
304  fhNCellsWithLargeTimeInCluster [il0][ism] = 0;
305  }
306  }
307 
308 // for(Int_t il0 = 0; il0 < 2; il0++)
309 // {
310 // for(Int_t i = 0; i < 7; i++)
311 // {
312 // fhEtaPhiLam0BinPtBinSMShared [il0][i] = 0 ;
313 // fhColRowLam0BinPtBinSMShared [il0][i] = 0 ;
314 // }
315 //
316 // for(Int_t ism =0; ism < 12; ism++)
317 // {
318 // fhTimeLam0BinPerSMShared[il0][ism] = 0;
319 // fhLam0PerSMShared [ism] = 0;
320 // fhLam1PerSMShared [ism] = 0;
321 // }
322 // }
323 
324  for(Int_t ism = 0; ism < 20; ism++)
325  {
326  fhLam0PerSM [ism] = 0;
327  fhLam1PerSM [ism] = 0;
330 // fhLam0PerSMSPDPileUp [ism] = 0;
331 // fhLam1PerSMSPDPileUp [ism] = 0;
332  }
333 
334  for(Int_t ilarge = 0; ilarge < 5; ilarge++)
335  {
338  }
339 
340  for(Int_t i = 0; i < 14; i++)
341  {
342  fhEBinClusterEtaPhi[i] = 0 ;
343  fhEBinClusterColRow[i] = 0 ;
344  fhEBinClusterEtaPhiPID[i] = 0 ;
345  fhEBinClusterColRowPID[i] = 0 ;
346  }
347 
348  for(Int_t igen = 0; igen < 10; igen++)
349  {
350  for(Int_t ip = 0; ip < fgkNGenTypes; ip++)
351  {
352  fhMergeGeneratorCluster [igen][ip] = 0;
353  fhMergeGeneratorClusterNotHijingBkg [igen][ip] = 0;
355  fhMergeGeneratorClusterHijingBkg [igen][ip] = 0;
356  fhCleanGeneratorCluster [igen][ip] = 0;
357 
363 
364  fhMergeGeneratorClusterEPrimRecoDiff [igen][ip] = 0;
368  fhCleanGeneratorClusterEPrimRecoDiff [igen][ip] = 0;
369  }
370  }
371 
372  for(Int_t icase = 0; icase < 6; icase++)
373  {
374  fhLocalRegionClusterEtaPhi [icase] = 0;
375  fhLocalRegionClusterEnergySum [icase] = 0;
387 
400  }
401 
402  // Initialize parameters
403  InitParameters();
404 }
405 
406 //_____________________________________________________________________
416 //__________________________________________________________________________
418  Float_t eta, Float_t phi,
419  Int_t mctag, TObjArray *clusterList)
420 {
421  Float_t radius = 0.2; // Hardcoded unless real use.
422 
423  // Accept cluster on EMCal limits - radius
424  // CAREFUL if used for Run2 with DCal.
425  if(phi < 3.15 - radius && phi > 1.4 + radius && TMath::Abs(eta) < 0.7-radius)
426  {
427  TString genName , genNameBkg ;
428  Int_t genIndex, genIndexBkg;
429  Int_t genBkgTag = -1;
430 
431  TString genName2 , genNameBkg2 ;
432  Int_t genIndex2, genIndexBkg2;
433  Int_t genBkgTag2 = -1;
434 
436  {
437  AliVCluster * calo = (AliVCluster*) (clusterList->At(icalo));
438 
439  genBkgTag = GetCocktailGeneratorBackgroundTag(calo, mctag, genName, genIndex, genNameBkg, genIndexBkg);
440 
441  }
442 
443  Bool_t mcDecay = kFALSE;
444  if( IsDataMC() && !genName.Contains("ijing") &&
447  mcDecay = kTRUE;
448 
449  fhLocalRegionClusterEtaPhi[0] ->Fill(eta,phi, GetEventWeight());
450 
451  Float_t sumE = 0;
452  Int_t sumM = 0;
453  Float_t sumEHi = 0;
454  Int_t sumMHi = 0;
455  Float_t sumEHi2 = 0;
456  Int_t sumMHi2 = 0;
457 
458  for(Int_t icalo2 = 0; icalo2 < clusterList->GetEntriesFast(); icalo2++)
459  {
460  if ( icalo2 == icalo ) continue;
461 
462  AliVCluster * calo2 = (AliVCluster*) (clusterList->At(icalo2));
463 
464  // Select clusters in a radius of R=0.2
465  calo2->GetMomentum(fMomentum2,GetVertex(0)) ;
466 
467  Float_t dEta = eta-fMomentum2.Eta();
468  Float_t dPhi = phi-GetPhi(fMomentum2.Phi());
469 
470  if(TMath::Abs(dPhi) >= TMath::Pi())
471  dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
472 
473  Float_t distance = TMath::Sqrt( dEta*dEta + dPhi*dPhi );
474 
475  genNameBkg2 = "";
476  genName2 = "";
477  genIndexBkg2 = -1;
478  genIndex2 = -1;
479  genBkgTag2 = -1;
480  if( IsDataMC() && IsStudyClusterOverlapsPerGeneratorOn() && calo2->GetNLabels() > 0 && distance < 0.4)
481  {
482  genBkgTag2 = GetCocktailGeneratorBackgroundTag(calo2, -1, genName2, genIndex2, genNameBkg2, genIndexBkg2);
483 
484  if( !genName.Contains("ijing") && !genName2.Contains("ijing")) fhDistance2AddedSignals ->Fill(en,distance,GetEventWeight());
485  if( !genName.Contains("ijing") && genName2.Contains("ijing")) fhDistanceAddedSignalsHijing->Fill(en,distance,GetEventWeight());
486  if( genName.Contains("ijing") && genName2.Contains("ijing")) fhDistance2Hijing ->Fill(en,distance,GetEventWeight());
487  }
488 
489  if ( distance > radius) continue;
490 
491  sumM++;
492  sumE += calo2->E();
493 
494  if( IsDataMC() && IsStudyClusterOverlapsPerGeneratorOn() && calo2->GetNLabels() > 0)
495  {
496  if(genName2.Contains("ijing"))
497  {
498  sumMHi++;
499  sumEHi += calo2->E();
500  }
501 
502  if(genName2.Contains("ijing") || genBkgTag2 == 1 || genBkgTag2 == 3)
503  {
504  sumMHi2++;
505  sumEHi2 += calo2->E();
506  }
507 
508  }
509  }
510 
511  Float_t sumMAdd = sumM-sumMHi;
512  Float_t sumEAdd = sumE-sumEHi;
513 
514  fhLocalRegionClusterEnergySum [0]->Fill(en,sumE,GetEventWeight());
516 
517  if(mcDecay)
518  {
521  }
522 
524  {
527 
528  if(mcDecay)
529  {
532  }
533  }
534 
536  {
537  fhLocalRegionClusterEnergySumHijing [0]->Fill(en,sumEHi ,GetEventWeight());
539 
542 
543  fhLocalRegionClusterEnergySumAdded [0]->Fill(en,sumEAdd,GetEventWeight());
545 
547  {
550 
553 
556  }
557 
558  if(mcDecay)
559  {
562 
565 
567  {
570 
573  }
574  }
575 
576  if(genBkgTag < 0) return;
577 
578  fhLocalRegionClusterEtaPhi [genBkgTag+1] ->Fill(eta,phi, GetEventWeight());
579 
580  fhLocalRegionClusterEnergySum [genBkgTag+1]->Fill(en,sumE,GetEventWeight());
581  fhLocalRegionClusterMultiplicity [genBkgTag+1]->Fill(en,sumM,GetEventWeight());
582 
583  fhLocalRegionClusterEnergySumHijing [genBkgTag+1]->Fill(en,sumEHi ,GetEventWeight());
584  fhLocalRegionClusterMultiplicityHijing[genBkgTag+1]->Fill(en,sumMHi ,GetEventWeight());
585 
586  fhLocalRegionClusterEnergySumAdded [genBkgTag+1]->Fill(en,sumEAdd,GetEventWeight());
587  fhLocalRegionClusterMultiplicityAdded [genBkgTag+1]->Fill(en,sumMAdd,GetEventWeight());
588 
590  {
593 
596 
599  }
600 
601  if(mcDecay)
602  {
603  fhLocalRegionClusterEnergySumMCPi0Decay [genBkgTag+1]->Fill(en,sumE,GetEventWeight());
604  fhLocalRegionClusterMultiplicityMCPi0Decay [genBkgTag+1]->Fill(en,sumM,GetEventWeight());
605 
606  fhLocalRegionClusterEnergySumHijingMCPi0Decay [genBkgTag+1]->Fill(en,sumEHi ,GetEventWeight());
607  fhLocalRegionClusterMultiplicityHijingMCPi0Decay[genBkgTag+1]->Fill(en,sumMHi ,GetEventWeight());
608 
609  fhLocalRegionClusterEnergySumAddedMCPi0Decay [genBkgTag+1]->Fill(en,sumEAdd,GetEventWeight());
610  fhLocalRegionClusterMultiplicityAddedMCPi0Decay [genBkgTag+1]->Fill(en,sumMAdd,GetEventWeight());
611 
613  {
616 
619 
622  }
623  }
624 
625  if(genBkgTag > 0) // Not clean, all merged
626  {
627  fhLocalRegionClusterEtaPhi [5] ->Fill(eta,phi, GetEventWeight());
628 
629  fhLocalRegionClusterEnergySum [5]->Fill(en,sumE,GetEventWeight());
631 
632  fhLocalRegionClusterEnergySumHijing [5]->Fill(en,sumEHi ,GetEventWeight());
634 
635  fhLocalRegionClusterEnergySumAdded [5]->Fill(en,sumEAdd,GetEventWeight());
637 
639  {
642 
645 
648  }
649 
650  if(mcDecay)
651  {
654 
657 
660 
662  {
665 
668 
671  }
672  }
673 
674  } // not clean, all merged
675 
676  } // check generator overlaps
677 
678  } // EMCal acceptance
679 }
680 
681 //_____________________________________________________________________
687 //__________________________________________________________________________
689 {
690  //
691  // Check the generators inside the cluster
692  TString genName = "", genNameBkg = "";
693  Int_t genIndex = -1, genIndexBkg = -1;
694  Int_t genBkgTag = GetCocktailGeneratorBackgroundTag(calo, mctag, genName, genIndex, genNameBkg, genIndexBkg);
695  if (genBkgTag == -1) return;
696  else if(genBkgTag > 3) printf("Bkg generator tag larger than 3\n");
697 
698  //
699  // Get primary particle info of main particle contributing to the cluster
700  Float_t eprim = 0;
701  //Float_t ptprim = 0;
702  Bool_t ok = kFALSE;
703  Int_t pdg = 0, status = 0, momLabel = -1;
704 
705  //fPrimaryMom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
706  fPrimaryMom = GetMCAnalysisUtils()->GetMother(calo->GetLabel(),GetMC(), pdg, status, ok, momLabel);
707 
708  if(ok)
709  {
710  eprim = fPrimaryMom.Energy();
711  //ptprim = fPrimaryMom.Pt();
712  }
713 
714  //
715  Float_t en = calo->E();
716 
717  if ( eprim < 0.1 || en < 0.5 ) return;
718 
719  //
720  Int_t partType = -1;
721  if ( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCPi0) ) partType = kmcGenPi0Merged;
722  else if( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCPi0Decay) ) partType = kmcGenPi0Decay;
723  else if( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCEtaDecay) ) partType = kmcGenEtaDecay;
724  else if( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCPhoton) ) partType = kmcGenPhoton;
725  else if( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCElectron) ) partType = kmcGenElectron;
726  else partType = kmcGenOther;
727 
728  Int_t genType = GetNCocktailGenNamesToCheck()-1;
729  for(Int_t igen = 1; igen < GetNCocktailGenNamesToCheck(); igen++)
730  {
731  if ( GetCocktailGenNameToCheck(igen).Contains(genName) &&
732  ( GetCocktailGenIndexToCheck(igen) < 0 || genIndex == GetCocktailGenIndexToCheck(igen)) )
733  {
734  genType = igen;
735  break;
736  }
737  }
738 
739  Float_t ratio = en / eprim;
740  Float_t diff = en - eprim;
741 
742  if ( genBkgTag > 0 )
743  {
744  fhMergeGeneratorCluster[0] [0] ->Fill(en,GetEventWeight());
745  fhMergeGeneratorCluster[0] [partType] ->Fill(en,GetEventWeight());
746  fhMergeGeneratorCluster[genType][0] ->Fill(en,GetEventWeight());
747  fhMergeGeneratorCluster[genType][partType] ->Fill(en,GetEventWeight());
748 
749  fhMergeGeneratorClusterEPrimRecoRatio[0] [0] ->Fill(en,ratio,GetEventWeight());
750  fhMergeGeneratorClusterEPrimRecoRatio[0] [partType] ->Fill(en,ratio,GetEventWeight());
751  fhMergeGeneratorClusterEPrimRecoRatio[genType][0] ->Fill(en,ratio,GetEventWeight());
752  fhMergeGeneratorClusterEPrimRecoRatio[genType][partType] ->Fill(en,ratio,GetEventWeight());
753 
754  fhMergeGeneratorClusterEPrimRecoDiff[0] [0] ->Fill(en,diff,GetEventWeight());
755  fhMergeGeneratorClusterEPrimRecoDiff[0] [partType] ->Fill(en,diff,GetEventWeight());
756  fhMergeGeneratorClusterEPrimRecoDiff[genType][0] ->Fill(en,diff,GetEventWeight());
757  fhMergeGeneratorClusterEPrimRecoDiff[genType][partType] ->Fill(en,diff,GetEventWeight());
758 
759  if(genBkgTag == 2)
760  {
762  fhMergeGeneratorClusterNotHijingBkg[0] [partType] ->Fill(en,GetEventWeight());
763  fhMergeGeneratorClusterNotHijingBkg[genType][0] ->Fill(en,GetEventWeight());
764  fhMergeGeneratorClusterNotHijingBkg[genType][partType] ->Fill(en,GetEventWeight());
765 
767  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[0] [partType] ->Fill(en,ratio,GetEventWeight());
769  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[genType][partType] ->Fill(en,ratio,GetEventWeight());
770 
774  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[genType][partType] ->Fill(en,diff,GetEventWeight());
775  }
776 
777  if(genBkgTag == 1)
778  {
780  fhMergeGeneratorClusterHijingBkg[0] [partType] ->Fill(en,GetEventWeight());
781  fhMergeGeneratorClusterHijingBkg[genType][0] ->Fill(en,GetEventWeight());
782  fhMergeGeneratorClusterHijingBkg[genType][partType] ->Fill(en,GetEventWeight());
783 
785  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[0] [partType] ->Fill(en,ratio,GetEventWeight());
786  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[genType][0] ->Fill(en,ratio,GetEventWeight());
787  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[genType][partType] ->Fill(en,ratio,GetEventWeight());
788 
790  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[0] [partType] ->Fill(en,diff,GetEventWeight());
792  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[genType][partType] ->Fill(en,diff,GetEventWeight());
793  }
794 
795  if ( genBkgTag == 3 )
796  {
798  fhMergeGeneratorClusterHijingAndOtherBkg[0] [partType] ->Fill(en,GetEventWeight());
800  fhMergeGeneratorClusterHijingAndOtherBkg[genType][partType] ->Fill(en,GetEventWeight());
801 
805  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[genType][partType] ->Fill(en,ratio,GetEventWeight());
806 
810  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[genType][partType] ->Fill(en,diff,GetEventWeight());
811  }
812  }
813  else
814  {
815  fhCleanGeneratorCluster[0] [0] ->Fill(en,GetEventWeight());
816  fhCleanGeneratorCluster[0] [partType] ->Fill(en,GetEventWeight());
817  fhCleanGeneratorCluster[genType][0] ->Fill(en,GetEventWeight());
818  fhCleanGeneratorCluster[genType][partType] ->Fill(en,GetEventWeight());
819 
820  fhCleanGeneratorClusterEPrimRecoRatio[0] [0] ->Fill(en,ratio,GetEventWeight());
821  fhCleanGeneratorClusterEPrimRecoRatio[0] [partType] ->Fill(en,ratio,GetEventWeight());
822  fhCleanGeneratorClusterEPrimRecoRatio[genType][0] ->Fill(en,ratio,GetEventWeight());
823  fhCleanGeneratorClusterEPrimRecoRatio[genType][partType] ->Fill(en,ratio,GetEventWeight());
824 
825  fhCleanGeneratorClusterEPrimRecoDiff[0] [0] ->Fill(en,diff,GetEventWeight());
826  fhCleanGeneratorClusterEPrimRecoDiff[0] [partType] ->Fill(en,diff,GetEventWeight());
827  fhCleanGeneratorClusterEPrimRecoDiff[genType][0] ->Fill(en,diff,GetEventWeight());
828  fhCleanGeneratorClusterEPrimRecoDiff[genType][partType] ->Fill(en,diff,GetEventWeight());
829  }
830 }
831 
832 
833 //_____________________________________________________________________
855 //_____________________________________________________________________
857  Int_t nMaxima, Bool_t matched, Int_t mctag, Float_t mcbin, Float_t egen, Int_t noverlaps)
858 {
859  Float_t ptcluster = fMomentum.Pt();
860  Float_t ecluster = fMomentum.E();
861  Float_t etacluster = fMomentum.Eta();
862  Float_t phicluster = fMomentum.Phi();
863 
864  if(phicluster < 0) phicluster+=TMath::TwoPi();
865 
866  AliDebug(2,Form("Current Event %d; Before selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f",
868  ecluster,ptcluster, phicluster*TMath::RadToDeg(),etacluster));
869 
870  fhClusterCutsE [1]->Fill( ecluster, GetEventWeight());
871  fhClusterCutsPt[1]->Fill(ptcluster, GetEventWeight());
872 
873  if(ecluster > 0.5) fhEtaPhi->Fill(etacluster, phicluster, GetEventWeight());
874 
875  if ( sm < fFirstModule || sm > fLastModule || sm >= fNModules || sm < 0 )
876  {
877  AliWarning(Form("SM %d out of expected range [%d, %d]",
878  sm, fFirstModule, fLastModule));
879  return kFALSE;
880  }
881 
882  fhEClusterSM ->Fill(ecluster , sm, GetEventWeight());
883  fhPtClusterSM->Fill(ptcluster, sm, GetEventWeight());
884 
885  //.......................................
886  //If too small or big energy, skip it
887  if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) return kFALSE ;
888 
889  AliDebug(2,Form("\t Cluster %d Pass E Cut",calo->GetID()));
890 
891  fhClusterCutsE [2]->Fill( ecluster, GetEventWeight());
892  fhClusterCutsPt[2]->Fill(ptcluster, GetEventWeight());
893 
894  //.......................................
895  // TOF cut, BE CAREFUL WITH THIS CUT
896  Double_t tof = calo->GetTOF()*1e9;
897  if(tof > 400) tof-=fConstantTimeShift; // only for MC, rest shift = 0
898 
899  if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
900 
901  AliDebug(2,Form("\t Cluster %d Pass Time Cut",calo->GetID()));
902 
903  fhClusterCutsE [3]->Fill( ecluster, GetEventWeight());
904  fhClusterCutsPt[3]->Fill(ptcluster, GetEventWeight());
905 
906  //.......................................
907  if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
908 
909  AliDebug(2,Form("\t Cluster %d Pass NCell Cut",calo->GetID()));
910 
911  fhClusterCutsE [4]->Fill( ecluster, GetEventWeight());
912  fhClusterCutsPt[4]->Fill(ptcluster, GetEventWeight());
913 
914  if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
915  AliDebug(2,Form("\t Cluster %d pass NLM %d of out of range",calo->GetID(), nMaxima));
916 
917  fhClusterCutsE [5]->Fill( ecluster, GetEventWeight());
918  fhClusterCutsPt[5]->Fill(ptcluster, GetEventWeight());
919 
920  //.......................................
921  //Check acceptance selection
922  if(IsFiducialCutOn())
923  {
925  if(! in ) return kFALSE ;
926  }
927 
928  AliDebug(2,Form("\t Fiducial cut passed"));
929 
930  fhClusterCutsE [6]->Fill( ecluster, GetEventWeight());
931  fhClusterCutsPt[6]->Fill(ptcluster, GetEventWeight());
932 
933  //.......................................
934  //Skip matched clusters with tracks
935 
936  if ( IsDataMC() )
937  {
938  // Fill selected cluster histograms in MC before TM
939  fhMCParticle[1]->Fill(ecluster, mcbin, GetEventWeight());
940 
941  if ( egen > 0.1 )
942  fhMCParticleVsErecEgen[1]->Fill(ecluster, mcbin, ecluster/egen, GetEventWeight());
943 
944  if ( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCConversion) )
945  fhMCParticleConverted[1]->Fill(ecluster, mcbin, GetEventWeight());
946 
947  fhMCParticleVsNOverlaps[1]->Fill(ecluster, mcbin, noverlaps, GetEventWeight());
948  }
949 
950  // Fill matching residual histograms before PID cuts
951  if ( fFillTMHisto )
952  FillTrackMatchingResidualHistograms(calo, 0, sm, matched, mctag, mcbin, egen, noverlaps);
953 
955  {
956  if(matched)
957  {
958  AliDebug(2,"\t Reject track-matched clusters");
959  return kFALSE ;
960  }
961  else
962  AliDebug(2,"\t Track-matching cut passed");
963  }// reject matched clusters
964 
965  fhClusterCutsE [7]->Fill( ecluster, GetEventWeight());
966  fhClusterCutsPt[7]->Fill(ptcluster, GetEventWeight());
967 
968  if ( IsDataMC() )
969  {
970  // Fill selected cluster histograms in MC before TM
971  fhMCParticle[2]->Fill(ecluster, mcbin, GetEventWeight());
972 
973  if ( egen > 0.1 )
974  fhMCParticleVsErecEgen[2]->Fill(ecluster, mcbin, ecluster/egen, GetEventWeight());
975 
976  if ( GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCConversion) )
977  fhMCParticleConverted[2]->Fill(ecluster, mcbin, GetEventWeight());
978 
979  fhMCParticleVsNOverlaps[2]->Fill(ecluster, mcbin, noverlaps, GetEventWeight());
980  }
981 
982  //.......................................
983  //Check Distance to Bad channel, set bit.
984  Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
985  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
986  if(distBad < fMinDist)
987  {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
988  return kFALSE ;
989  }
990  else AliDebug(2,Form("\t Bad channel cut passed %4.2f > %2.2f",distBad, fMinDist));
991 
992  fhClusterCutsE [8]->Fill( ecluster, GetEventWeight());
993  fhClusterCutsPt[8]->Fill(ptcluster, GetEventWeight());
994 
995  AliDebug(1,Form("Current Event %d; After selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f",
997  ecluster, ptcluster,fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
998 
999  //All checks passed, cluster selected
1000  return kTRUE;
1001 
1002 }
1003 
1004 //___________________________________________
1009 //___________________________________________
1011 {
1012  if( !GetMC() )
1013  {
1014  AliFatal("MCEvent not available, is the MC handler called? STOP");
1015  return;
1016  }
1017 
1018  Double_t photonY = -100 ;
1019  Double_t photonE = -1 ;
1020  Double_t photonPt = -1 ;
1021  Double_t photonPhi = 100 ;
1022  Double_t photonEta = -1 ;
1023 
1024  Int_t pdg = 0 ;
1025  Int_t tag = 0 ;
1026  Int_t status = 0 ;
1027  Int_t mcIndex = 0 ;
1028  Int_t nprim = GetMC()->GetNumberOfTracks();
1029  Bool_t inacceptance = kFALSE ;
1030 
1031  AliVParticle * primary = 0;
1032 
1033  for(Int_t i=0 ; i < nprim; i++)
1034  {
1035  if ( !GetReader()->AcceptParticleMCLabel( i ) ) continue ;
1036 
1037  primary = GetMC()->GetTrack(i);
1038  if(!primary)
1039  {
1040  AliWarning("primaries pointer not available!!");
1041  continue;
1042  }
1043 
1044  pdg = primary->PdgCode();
1045  status = primary->MCStatusCode();
1046 
1047  // Protection against floating point exception
1048  if ( primary->E() == TMath::Abs(primary->Pz()) ||
1049  (primary->E() - primary->Pz()) < 1e-3 ||
1050  (primary->E() + primary->Pz()) < 0 ) continue ;
1051 
1052  //printf("i %d, %s %d %s %d \n",i, GetMC()->Particle(i)->GetName(), GetMC()->Particle(i)->GetPdgCode(),
1053  // prim->GetName(), prim->GetPdgCode());
1054 
1055  // Photon kinematics
1056  primary->Momentum(fMomentum);
1057 
1058  photonY = 0.5*TMath::Log((primary->E()+primary->Pz())/(primary->E()-primary->Pz())) ;
1059 
1060  //printf("particle %d, pdg %d, status %d, phys prim %d, pT %2.2f\n",i,pdg,status,primary->IsPhysicalPrimary(),primary->Pt());
1061 
1062  // Select only photons in the final state
1063  if(pdg != 22 ) continue ;
1064 
1065  // If too small or too large pt, skip, same cut as for data analysis
1066  photonPt = fMomentum.Pt () ;
1067 
1068  if(photonPt < GetMinPt() || photonPt > GetMaxPt() ) continue ;
1069 
1070  photonE = fMomentum.E () ;
1071  photonEta = fMomentum.Eta() ;
1072  photonPhi = fMomentum.Phi() ;
1073 
1074  if(photonPhi < 0) photonPhi+=TMath::TwoPi();
1075 
1076  // Check if photons hit desired acceptance
1077  inacceptance = kTRUE;
1078 
1079  // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
1080  if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter())) inacceptance = kFALSE ;
1081 
1082  // Check if photons hit the Calorimeter acceptance
1083  if(IsRealCaloAcceptanceOn()) // defined on base class
1084  {
1085  if(!GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(GetCalorimeter(), primary)) inacceptance = kFALSE ;
1086  }
1087 
1088  // Get tag of this particle photon from fragmentation, decay, prompt ...
1089  // Set the origin of the photon.
1090  tag = GetMCAnalysisUtils()->CheckOrigin(i, GetMC(),
1091  GetReader()->GetNameOfMCEventHederGeneratorToAccept(),
1092  photonE); // not needed
1093 
1094  if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
1095  {
1096  // A conversion photon from a hadron, skip this kind of photon
1097  // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
1098  // GetMCAnalysisUtils()->PrintMCTag(tag);
1099 
1100  continue;
1101  }
1102 
1103  //if(fStudyActivityNearCluster && photonE > 5) // cut at 5 GeV to avoid too many iterations
1104  // DistanceToAddedSignalAtGeneratorLevel(i,nprim, GetMC(), photonE,photonEta,photonPhi);
1105 
1106  // Consider only final state particles, but this depends on generator,
1107  // status 1 is the usual one, in case of not being ok, leave the possibility
1108  // to not consider this.
1109  if(status > 1) continue ; // Avoid "partonic" photons
1110 
1111  Bool_t takeIt = kFALSE ;
1112  if(status == 1 && GetMCAnalysisUtils()->GetMCGenerator() != AliMCAnalysisUtils::kBoxLike ) takeIt = kTRUE ;
1113 
1114  if (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) continue;
1115 
1116  // Origin of photon
1117  if (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
1118  {
1119  mcIndex = kmcPPrompt;
1120  }
1121  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
1122  {
1123  mcIndex = kmcPFragmentation ;
1124  }
1125  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
1126  {
1127  mcIndex = kmcPISR;
1128  }
1129  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
1130  {
1131  mcIndex = kmcPPi0Decay;
1132  }
1133  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
1134  {
1135  mcIndex = kmcPEtaDecay;
1136  }
1137  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
1138  {
1139  mcIndex = kmcPOtherDecay;
1140  }
1141  else
1142  {
1143  // Other decay but from non final state particle
1144  mcIndex = kmcPOtherDecay;
1145  } // Other origin
1146 
1147  if(!takeIt && (mcIndex == kmcPPi0Decay || mcIndex == kmcPOtherDecay)) takeIt = kTRUE ;
1148 
1149  if(!takeIt) continue ;
1150 
1151  // Fill histograms for all photons
1152  fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY, GetEventWeight()) ;
1153  if(TMath::Abs(photonY) < 1.0)
1154  {
1155  fhEPrimMC [kmcPPhoton]->Fill(photonE , GetEventWeight()) ;
1156  fhPtPrimMC [kmcPPhoton]->Fill(photonPt, GetEventWeight()) ;
1157  fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi, GetEventWeight()) ;
1158  fhEtaPrimMC[kmcPPhoton]->Fill(photonE , photonEta, GetEventWeight()) ;
1159  }
1160 
1161  if(inacceptance)
1162  {
1163  fhEPrimMCAcc [kmcPPhoton]->Fill(photonE , GetEventWeight()) ;
1164  fhPtPrimMCAcc [kmcPPhoton]->Fill(photonPt, GetEventWeight()) ;
1165  fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi, GetEventWeight()) ;
1166  fhEtaPrimMCAcc[kmcPPhoton]->Fill(photonE , photonEta, GetEventWeight()) ;
1167  fhYPrimMCAcc [kmcPPhoton]->Fill(photonE , photonY , GetEventWeight()) ;
1168  } // Accepted
1169 
1170  // Fill histograms for photons origin
1171  if(mcIndex < fNPrimaryHistograms)
1172  {
1173  fhYPrimMC[mcIndex]->Fill(photonPt, photonY, GetEventWeight()) ;
1174  if(TMath::Abs(photonY) < 1.0)
1175  {
1176  fhEPrimMC [mcIndex]->Fill(photonE , GetEventWeight()) ;
1177  fhPtPrimMC [mcIndex]->Fill(photonPt, GetEventWeight()) ;
1178  fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi, GetEventWeight()) ;
1179  fhEtaPrimMC[mcIndex]->Fill(photonE , photonEta, GetEventWeight()) ;
1180  }
1181 
1182  if(inacceptance)
1183  {
1184  fhEPrimMCAcc [mcIndex]->Fill(photonE , GetEventWeight()) ;
1185  fhPtPrimMCAcc [mcIndex]->Fill(photonPt, GetEventWeight()) ;
1186  fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi, GetEventWeight()) ;
1187  fhEtaPrimMCAcc[mcIndex]->Fill(photonE , photonEta, GetEventWeight()) ;
1188  fhYPrimMCAcc [mcIndex]->Fill(photonE , photonY , GetEventWeight()) ;
1189  } // Accepted
1190  }
1191  } // loop on primaries
1192 }
1193 
1194 //________________________________________________________________________________
1199 //________________________________________________________________________________
1200 //void AliAnaPhoton::DistanceToAddedSignal(Int_t label, Int_t nprim,
1201 // Float_t photonE, Float_t photonEta, Float_t photonPhi)
1202 //{
1203 // //Double_t addedY = -100 ;
1204 // //Double_t addedE = -1 ;
1205 // Double_t addedPt = -1 ;
1206 // Double_t addedPhi = 100 ;
1207 // Double_t addedEta = -1 ;
1208 //
1209 // //Int_t pdg = 0 ;
1210 // Int_t status = 0 ;
1211 // Bool_t inacceptance = kFALSE ;
1212 //
1213 // Int_t pdgMomOrg = 0, statusMomOrg = -1, momLabelOrg = -1;
1214 // Bool_t okMomOrg = kFALSE;
1215 // GetMCAnalysisUtils()->GetMother(label,GetMC(), pdgMomOrg, statusMomOrg, okMomOrg, momLabelOrg);
1216 //
1217 // AliVParticle * primary = 0;
1218 // AliVParticle * primaryMom = 0;
1219 //
1220 // TString genName;
1221 // (GetMC())->GetCocktailGenerator(label,genName);
1222 //
1223 // for(Int_t i=0 ; i < nprim; i++)
1224 // {
1225 // if( i == label ) continue ;
1226 //
1227 // if ( !GetReader()->AcceptParticleMCLabel( i ) ) continue ;
1228 //
1229 // //printf("particle %d\n",i);
1230 //
1231 // Int_t pdgMom = 0, statusMom = -1, momLabel = -1;
1232 // Bool_t okMom = kFALSE;
1233 // GetMCAnalysisUtils()->GetMother(i,GetMC(), pdgMom, statusMom, okMom, momLabel);
1234 // //printf("\t mom label %d\n",momLabel);
1235 //
1236 // if(momLabel < 0) continue;
1237 //
1238 // Int_t grandMomLabel = -1;
1239 // primary = GetMC()->Particle(i) ;
1240 // if(!primary)
1241 // {
1242 // AliWarning("primaries pointer not available!!");
1243 // continue;
1244 // }
1245 //
1246 // //pdg = primary->GetPdgCode();
1247 // status = primary->GetStatusCode();
1248 //
1249 // // Protection against floating point exception
1250 // if ( primary->E() == TMath::Abs(primary->Pz()) ||
1251 // (primary->E() - primary->Pz()) < 1e-3 ||
1252 // (primary->E() + primary->Pz()) < 0 ) continue ;
1253 //
1254 // //printf("i %d, %s %d %s %d \n",i, GetMC()->Particle(i)->GetName(), GetMC()->Particle(i)->GetPdgCode(),
1255 // // prim->GetName(), prim->GetPdgCode());
1256 //
1257 // // Photon kinematics
1258 // primary->Momentum(fMomentum);
1259 //
1260 // //addedY = 0.5*TMath::Log((primary->E+primary->Pz())/(primary->E-primary->Pz())) ;
1261 //
1262 // if(momLabel >= 0)
1263 // {
1264 // primaryMom = GetMC()->Particle(momLabel) ;
1265 // grandMomLabel = primaryMom->GetFirstMother();
1266 // }
1267 // //printf("\t grandmom %d\n",grandMomLabel);
1268 //
1269 // // Select only photons in the final state
1270 // //if(pdg != 22 ) continue ;
1271 //
1272 // // If too small or too large pt, skip, same cut as for data analysis
1273 // addedPt = fMomentum.Pt () ;
1274 //
1275 // if(addedPt < GetMinPt() || addedPt > GetMaxPt() ) continue ;
1276 //
1277 // //addedE = fMomentum.E () ;
1278 // addedEta = fMomentum.Eta() ;
1279 // addedPhi = fMomentum.Phi() ;
1280 //
1281 // if(addedPhi < 0) addedPhi+=TMath::TwoPi();
1282 //
1283 // // Check if photons hit desired acceptance
1284 // inacceptance = kTRUE;
1285 //
1286 // // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
1287 // if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter())) inacceptance = kFALSE ;
1288 //
1289 // // Check if photons hit the Calorimeter acceptance
1290 // if(IsRealCaloAcceptanceOn()) // defined on base class
1291 // {
1292 // if(!GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(GetCalorimeter(), primary)) inacceptance = kFALSE ;
1293 // }
1294 // //printf("\t in acceptance %d", inacceptance);
1295 // if(!inacceptance) continue;
1296 //
1297 // TString genName2;
1298 // (GetMC())->GetCocktailGenerator(i,genName2);
1299 //
1300 // Float_t dEta = photonEta-addedEta;
1301 // Float_t dPhi = photonPhi-addedPhi;
1302 // Float_t distance = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1303 //
1304 // if (( momLabel != momLabelOrg && momLabelOrg >= 0) || momLabelOrg < 0)
1305 // {
1306 // if(!genName2.Contains("ijing"))
1307 // {
1308 // if(grandMomLabel < 0)
1309 // {
1310 // if ( !genName.Contains("ijing") ) fhDistanceAddedPhotonAddedPrimarySignal ->Fill(photonE,distance,GetEventWeight());
1311 // if ( genName.Contains("ijing") ) fhDistanceHijingPhotonAddedPrimarySignal ->Fill(photonE,distance,GetEventWeight());
1312 // }
1313 // else
1314 // {
1315 // if ( !genName.Contains("ijing") ) fhDistanceAddedPhotonAddedSecondarySignal ->Fill(photonE,distance,GetEventWeight());
1316 // if ( genName.Contains("ijing") ) fhDistanceHijingPhotonAddedSecondarySignal->Fill(photonE,distance,GetEventWeight());
1317 // }
1318 // }
1319 // }
1320 //
1321 // if(!genName.Contains("ijing") && genName2.Contains("ijing") && status == 0) fhDistanceAddedPhotonHijingSecondary->Fill(photonE,distance,GetEventWeight());
1322 // }
1323 //}
1324 //
1325 //________________________________________________________________________________
1327 //________________________________________________________________________________
1328 void AliAnaPhoton::FillPileUpHistograms(AliVCluster* cluster, AliVCaloCells *cells,
1329  Int_t absIdMax)
1330 {
1331  Float_t pt = fMomentum.Pt();
1332  Float_t time = cluster->GetTOF()*1.e9;
1333  if ( time > 400 ) time-=fConstantTimeShift; // Only for MC, rest shift = 0
1334 
1335  AliVEvent * event = GetReader()->GetInputEvent();
1336 
1337  if(GetReader()->IsPileUpFromSPD()) fhPtPhotonPileUp[0]->Fill(pt, GetEventWeight()) ;
1338  if(GetReader()->IsPileUpFromEMCal()) fhPtPhotonPileUp[1]->Fill(pt, GetEventWeight()) ;
1339  if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPhotonPileUp[2]->Fill(pt, GetEventWeight()) ;
1340  if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPhotonPileUp[3]->Fill(pt, GetEventWeight()) ;
1341  if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPhotonPileUp[4]->Fill(pt, GetEventWeight()) ;
1342  if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPhotonPileUp[5]->Fill(pt, GetEventWeight()) ;
1343  if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPhotonPileUp[6]->Fill(pt, GetEventWeight()) ;
1344 
1345  fhTimePtPhotonNoCut->Fill(pt, time, GetEventWeight());
1346  if(GetReader()->IsPileUpFromSPD()) fhTimePtPhotonSPD->Fill(pt, time, GetEventWeight());
1347 
1348  //
1349  // Cells inside the cluster
1350  //
1351 
1352  // Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
1353  if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(time) < 30)
1354  {
1355  for (Int_t ipos = 0; ipos < cluster->GetNCells(); ipos++)
1356  {
1357  Int_t absId = cluster->GetCellsAbsId()[ipos];
1358 
1359  if( absId == absIdMax ) continue ;
1360 
1361  Double_t tcell = cells->GetCellTime(absId);
1362  Float_t amp = cells->GetCellAmplitude(absId);
1363  Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
1364 
1365  GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,tcell,cells);
1366  tcell*=1e9;
1367  if(tcell > 400)tcell-=fConstantTimeShift;
1368 
1369  Float_t diff = (time-tcell);
1370 
1371  if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
1372 
1380 
1381  } // loop
1382  }
1383 
1384  AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
1385  AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
1386 
1387  //
1388  // N pile up vertices
1389  //
1390 
1391  Int_t nVtxSPD = -1;
1392  Int_t nVtxTrk = -1;
1393 
1394  if (esdEv)
1395  {
1396  nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
1397  nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
1398 
1399  } // ESD
1400  else if (aodEv)
1401  {
1402  nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
1403  nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
1404  } // AOD
1405 
1406  if(pt < 8)
1407  {
1408  fhTimeNPileUpVertSPD ->Fill(time, nVtxSPD, GetEventWeight());
1409  fhTimeNPileUpVertTrack->Fill(time, nVtxTrk, GetEventWeight());
1410  }
1411 
1412  fhPtPhotonNPileUpSPDVtx->Fill(pt, nVtxSPD, GetEventWeight());
1413  fhPtPhotonNPileUpTrkVtx->Fill(pt, nVtxTrk, GetEventWeight());
1414 
1415  if(TMath::Abs(time) < 25)
1416  {
1417  fhPtPhotonNPileUpSPDVtxTimeCut->Fill(pt, nVtxSPD, GetEventWeight());
1418  fhPtPhotonNPileUpTrkVtxTimeCut->Fill(pt, nVtxTrk, GetEventWeight());
1419  }
1420 
1421  if(time < 75 && time > -25)
1422  {
1423  fhPtPhotonNPileUpSPDVtxTimeCut2->Fill(pt, nVtxSPD, GetEventWeight());
1424  fhPtPhotonNPileUpTrkVtxTimeCut2->Fill(pt, nVtxTrk, GetEventWeight());
1425  }
1426 }
1427 
1428 //_________________________________________________________________________________
1430 //_________________________________________________________________________________
1431 void AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, Int_t sm, Int_t mcTag,
1432  Int_t nlm, Bool_t matched,
1433  Float_t maxCellFraction, Int_t & largeTime)
1434 {
1435  if(!fFillSSHistograms || GetMixedEvent()) return;
1436 
1437  Float_t energy = cluster->E();
1438  Int_t ncells = cluster->GetNCells();
1439  Float_t lambda0 = cluster->GetM02();
1440  Float_t lambda1 = cluster->GetM20();
1441  Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
1442 
1443  Float_t pt = fMomentum.Pt();
1444  Float_t eta = fMomentum.Eta();
1445  Float_t phi = GetPhi(fMomentum.Phi());
1446 
1447  fhLam0E ->Fill(energy, lambda0, GetEventWeight());
1448  fhLam0Pt->Fill(pt , lambda0, GetEventWeight());
1449  fhLam1E ->Fill(energy, lambda1, GetEventWeight());
1450  fhLam1Pt->Fill(pt , lambda1, GetEventWeight());
1451 
1453  {
1454  fhLam0PerSM[sm]->Fill(pt, lambda0, GetEventWeight());
1455  fhLam1PerSM[sm]->Fill(pt, lambda1, GetEventWeight());
1456  }
1457 
1459  {
1460  fhDispE ->Fill(energy, disp , GetEventWeight());
1461  fhDispPt->Fill(pt , disp , GetEventWeight());
1462  }
1463 
1465  {
1466  if(nlm==1)
1467  {
1468  fhLam0PtNLM1->Fill(pt, lambda0, GetEventWeight());
1469  fhLam1PtNLM1->Fill(pt, lambda1, GetEventWeight());
1470  }
1471  else if(nlm==2)
1472  {
1473  fhLam0PtNLM2->Fill(pt, lambda0, GetEventWeight());
1474  fhLam1PtNLM2->Fill(pt, lambda1, GetEventWeight());
1475  }
1476  }
1477 
1478  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
1479  sm >= GetFirstSMCoveredByTRD() )
1480  {
1481  fhLam0ETRD ->Fill(energy, lambda0, GetEventWeight());
1482  fhLam0PtTRD->Fill(pt , lambda0, GetEventWeight());
1483  fhLam1ETRD ->Fill(energy, lambda1, GetEventWeight());
1485  fhDispETRD ->Fill(energy, disp, GetEventWeight());
1486  }
1487 
1488  //
1489  // EMCal SM regions
1490  //
1491  if(cluster->IsEMCAL() && fFillEMCALRegionSSHistograms)
1492  {
1493 // Bool_t shared = GetCaloUtils()->IsClusterSharedByTwoSuperModules(GetEMCALGeometry(),cluster);
1494 // if(GetReader()->IsPileUpFromSPD())
1495 // {
1496 // fhLam0PerSMSPDPileUp[sm]->Fill(pt, lambda0, GetEventWeight());
1497 // fhLam1PerSMSPDPileUp[sm]->Fill(pt, lambda1, GetEventWeight());
1498 // }
1499 
1500  Int_t etaRegion = -1, phiRegion = -1;
1501  GetCaloUtils()->GetEMCALSubregion(cluster,GetReader()->GetEMCALCells(),etaRegion,phiRegion);
1502  if(etaRegion >= 0 && etaRegion < 4 && phiRegion >=0 && phiRegion < 3)
1503  {
1504 
1505  // fhLam0EMCALRegion[etaRegion][phiRegion]->Fill(pt,lambda0, GetEventWeight());
1506  //
1507  // if(GetFirstSMCoveredByTRD() >= 0 && sm >= GetFirstSMCoveredByTRD() )
1508  // {
1509  // fhLam0EMCALRegionTRD[etaRegion][phiRegion]->Fill(pt, lambda0, GetEventWeight());
1510  // }
1511 
1512  fhLam0EMCALRegionPerSM[etaRegion][phiRegion][sm]->Fill(pt,lambda0, GetEventWeight());
1513  fhLam1EMCALRegionPerSM[etaRegion][phiRegion][sm]->Fill(pt,lambda1, GetEventWeight());
1514 
1515 // if(shared)
1516 // {
1517 // fhLam0PerSMShared[sm]->Fill(pt,lambda0, GetEventWeight());
1518 // fhLam1PerSMShared[sm]->Fill(pt,lambda1, GetEventWeight());
1519 // }
1520  }
1521 
1522  // Define l0 bin and cluster pt bin
1523  Int_t l0bin = -1;
1524  if ( lambda0 >=0.30 && lambda0 <= 0.40 ) l0bin = 1;
1525  else if( lambda0 >=0.23 && lambda0 <= 0.26 ) l0bin = 0;
1526 
1527  Float_t ptLimit[] = {2,3,4,5,6,8,10,12};
1528  Int_t ptbin = -1;
1529  for(Int_t ipt = 0; ipt < 7; ipt++)
1530  {
1531  if( pt >= ptLimit[ipt] && pt < ptLimit[ipt+1] )
1532  {
1533  ptbin = ipt;
1534  break;
1535  }
1536  }
1537 
1538  if(l0bin!=-1)
1539  {
1540  fhLam1Lam0BinPerSM[l0bin][sm]->Fill(pt,lambda1);
1541 
1542  if(ptbin >= 0)
1543  {
1544  fhEtaPhiLam0BinPtBin[l0bin][ptbin]->Fill(eta, phi, GetEventWeight());
1545 // if(shared) fhEtaPhiLam0BinPtBinSMShared[l0bin][ptbin]->Fill(eta, phi, GetEventWeight());
1546  }
1547  }
1548 
1549  //
1550  // Histograms with time of secondary cells, only if they contribute with a non null
1551  // weight to the shower shape
1552  //
1553  AliVCaloCells* cells = GetReader()->GetEMCALCells();
1554 
1555  Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
1556 
1557  //
1558  // Sort the cluster cells in energy, find maximum energy cell
1559  //
1560  Int_t ncell = cluster->GetNCells();
1561  Int_t * sortList = new Int_t [ncell];
1562  Float_t * enerList = new Float_t [ncell];
1563  Double_t * timeList = new Double_t[ncell];
1564 
1565  for(Int_t icell = 0; icell < ncell; icell++)
1566  {
1567  Int_t absId = cluster->GetCellAbsId(icell);
1568  Float_t cellE = cells->GetCellAmplitude(absId);
1569  Double_t cellTime = cells->GetCellTime(absId);
1570 
1571  GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,cellE,cellTime,cells);
1572  cellTime*=1e9;
1573  if(cellTime > 400)cellTime-=fConstantTimeShift;
1574 
1575  enerList[icell] = cellE;
1576  timeList[icell] = cellTime;
1577  }
1578 
1579  TMath::Sort(ncell,enerList,sortList);
1580 
1581  //
1582  // Location of max cell
1583  // First in sorted list.
1584  //
1585  Int_t absIdMax = cluster->GetCellAbsId(sortList[0]);
1586  // Float_t cellEMax = enerList[sortList[0]];
1587  Double_t cellTimeMax = timeList[sortList[0]];
1588 
1589  //
1590  // Cluster cell loop, select only secondary cells with enough contribution to cluster
1591  // Start from second highest energy cell
1592  //
1593  largeTime = 0;
1594  Float_t largeTimeE = 0;
1595  for(Int_t icell = 1; icell < ncell; icell++)
1596  {
1597  Int_t absId = cluster->GetCellAbsId(sortList[icell]);
1598  Float_t cellE = enerList[sortList[icell]];
1599  Double_t cellTime = timeList[sortList[icell]];
1600 
1601  //printf("\t icell %d, absId %d, ener %2.2f, time %2.2f, sort %d\n",icell,absId,cellE,cellTime,sortList[icell]);
1602 
1603  if ( absId == absIdMax ) continue;
1604 
1605  Float_t weight = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cellE,cluster->E());
1606 
1607  //printf("Cluster E %2.2f, cell E %2.2f, time %2.2f, weight %2.3f\n",cluster->E(),cellE,cellTime,weight);
1608  if(weight < 0.01) continue;
1609 
1610  if ( TMath::Abs(cellTime) > 50 )
1611  {
1612  largeTime++;
1613  largeTimeE += cellE;
1614  //printf(" cluster pT %2.2f, cell E %2.2f, large time %d; ",pt, cellE, largeTime);
1615  }
1616  if ( l0bin == -1 ) continue;
1617 
1618  fhTimeLam0BinPerSM [l0bin][sm]->Fill(pt, cellTime, GetEventWeight());
1619  fhTimeLam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellTime, weight*GetEventWeight());
1620 
1621  fhDTimeLam0BinPerSM [l0bin][sm]->Fill(pt, cellTimeMax-cellTime, GetEventWeight());
1622  fhDTimeLam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellTimeMax-cellTime, weight*GetEventWeight());
1623 
1624  fhCellClusterEFracLam0BinPerSM [l0bin][sm]->Fill(pt, cellE/cluster->E(), GetEventWeight());
1625 // fhCellClusterEFracLam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellE/cluster->E(), weight*GetEventWeight());
1626 
1627  fhCellClusterELam0BinPerSM [l0bin][sm]->Fill(pt, cellE, GetEventWeight());
1628  fhCellClusterELam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellE, weight*GetEventWeight());
1629 
1630  if(ptbin >= 0)
1631  {
1632  fhCellClusterIndexEAndTime[l0bin][ptbin] ->Fill(cellTime, icell , GetEventWeight());
1633  fhCellClusterEAndTime [l0bin][ptbin] ->Fill(cellTime, cellE , GetEventWeight());
1634  fhCellClusterEFracAndTime [l0bin][ptbin] ->Fill(cellTime, cellE/cluster->E(), GetEventWeight());
1635  }
1636 
1637  if ( TMath::Abs(cellTime) > 50 )
1638  {
1639  fhCellClusterEFracLam0BinPerSMLargeTime [l0bin][sm]->Fill(pt, cellE/cluster->E(), GetEventWeight());
1640  fhCellClusterELam0BinPerSMLargeTime [l0bin][sm]->Fill(pt, cellE , GetEventWeight());
1641  fhCellClusterIndexELam0BinPerSMLargeTime[l0bin][sm]->Fill(pt, icell , GetEventWeight());
1642  }
1643 // if(shared)
1644 // fhTimeLam0BinPerSMShared[l0bin][sm]->Fill(pt, cellTime, GetEventWeight());
1645 
1646  if(ptbin >= 0)
1647  {
1648  Int_t icol = -1, icolAbs = -1;
1649  Int_t irow = -1, irowAbs = -1;
1650  Int_t iRCU = -1;
1651  GetModuleNumberCellIndexesAbsCaloMap(absId,GetCalorimeter(), icol, irow, iRCU, icolAbs, irowAbs);
1652 
1653  fhColRowLam0BinPtBin [l0bin][ptbin]->Fill(icolAbs, irowAbs, GetEventWeight());
1654  fhColRowLam0BinPtBinWeighted[l0bin][ptbin]->Fill(icolAbs, irowAbs, weight*GetEventWeight());
1655 
1656  if ( TMath::Abs(cellTime) > 50 )
1657  fhColRowLam0BinPtBinLargeTime[l0bin][ptbin]->Fill(icolAbs, irowAbs, GetEventWeight());
1658 
1659 // if(shared) fhColRowLam0BinPtBinSMShared[l0bin][ptbin]->Fill(icolAbs, irowAbs, GetEventWeight());
1660  }
1661 
1662 
1663  } // cluster cell loop
1664  //if(largeTime > 0) printf("\n");
1665 
1666  if ( largeTime > 0 )
1667  {
1668  if ( l0bin != -1 )
1669  {
1670  if ( ptbin > 0 )
1671  fhEtaPhiLargeTimeInClusterCell[l0bin][ptbin]->Fill(eta, phi, GetEventWeight());
1672 
1673  fhCellClusterEFracLam0BinPerSMLargeTimeTotal[l0bin][sm]->Fill(pt, largeTimeE/cluster->E(), GetEventWeight());
1674  fhNCellsWithLargeTimeInCluster [l0bin][sm]->Fill(pt, largeTime , GetEventWeight());
1675  }
1676 
1677  fhLam0PerSMLargeTimeInClusterCell[sm]->Fill(pt, lambda0, GetEventWeight());
1678  fhLam1PerSMLargeTimeInClusterCell[sm]->Fill(pt, lambda1, GetEventWeight());
1679 
1680  if(largeTime < 6)
1681  {
1682  fhLam0PerNLargeTimeInClusterCell[largeTime-1]->Fill(pt, lambda0, GetEventWeight());
1683  fhLam1PerNLargeTimeInClusterCell[largeTime-1]->Fill(pt, lambda1, GetEventWeight());
1684  }
1685  }
1686 
1687  delete [] sortList ;
1688  delete [] enerList ;
1689  delete [] timeList ;
1690  //printf("Cluster %d, E %2.2f, sm %d, eta %2.2f, phi %2.2f ---> region %d %d\n",cluster->GetID(),cluster->E(),sm,eta,RadToDeg(phi),etaRegion,phiRegion);
1691  } // Check EMCal SM and regions
1692 
1693  Float_t l0 = 0., l1 = 0.;
1694  Float_t dispp= 0., dEta = 0., dPhi = 0.;
1695  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1697  {
1699  l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
1700  //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",
1701  // l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
1702  //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
1703  // disp, dPhi+dEta );
1704  fhDispEtaE -> Fill(energy, dEta , GetEventWeight());
1705  fhDispPhiE -> Fill(energy, dPhi , GetEventWeight());
1706  fhSumEtaE -> Fill(energy, sEta , GetEventWeight());
1707  fhSumPhiE -> Fill(energy, sPhi , GetEventWeight());
1708  fhSumEtaPhiE -> Fill(energy, sEtaPhi , GetEventWeight());
1709  fhDispEtaPhiDiffE -> Fill(energy, dPhi-dEta, GetEventWeight());
1710 
1711  if(dEta+dPhi>0)fhSphericityE -> Fill(energy,(dPhi-dEta)/(dEta+dPhi) , GetEventWeight());
1712  if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.), GetEventWeight());
1713  if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.), GetEventWeight());
1714 
1715  Int_t ebin = -1;
1716  if (energy < 2 ) ebin = 0;
1717  else if (energy < 4 ) ebin = 1;
1718  else if (energy < 6 ) ebin = 2;
1719  else if (energy < 10) ebin = 3;
1720  else if (energy < 15) ebin = 4;
1721  else if (energy < 20) ebin = 5;
1722  else ebin = 6;
1723 
1724  fhDispEtaDispPhi[ebin]->Fill(dEta , dPhi, GetEventWeight());
1725  fhLambda0DispEta[ebin]->Fill(lambda0, dEta, GetEventWeight());
1726  fhLambda0DispPhi[ebin]->Fill(lambda0, dPhi, GetEventWeight());
1727  }
1728 
1729  // If track-matching was off, check effect of track-matching cut
1730  //
1731  if ( !fRejectTrackMatch && matched )
1732  {
1733  fhLam0ETM ->Fill(energy, lambda0, GetEventWeight());
1734  fhLam0PtTM->Fill(pt , lambda0, GetEventWeight());
1735  fhLam1ETM ->Fill(energy, lambda1, GetEventWeight());
1737  fhDispETM ->Fill(energy, disp , GetEventWeight());
1738 
1739  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
1740  sm >= GetFirstSMCoveredByTRD() )
1741  {
1742  fhLam0ETMTRD ->Fill(energy, lambda0, GetEventWeight());
1743  fhLam0PtTMTRD->Fill(pt , lambda0, GetEventWeight());
1744  fhLam1ETMTRD ->Fill(energy, lambda1, GetEventWeight());
1746  fhDispETMTRD ->Fill(energy, disp , GetEventWeight());
1747  }
1748 
1749  } // If track-matching was off, check effect of matching residual cut
1750 
1752  {
1753  if(energy < 2)
1754  {
1755  fhNCellsLam0LowE ->Fill(ncells, lambda0, GetEventWeight());
1756  fhNCellsLam1LowE ->Fill(ncells, lambda1, GetEventWeight());
1757  fhNCellsDispLowE ->Fill(ncells, disp , GetEventWeight());
1758 
1759  fhLam1Lam0LowE ->Fill(lambda1, lambda0, GetEventWeight());
1760  fhLam0DispLowE ->Fill(lambda0, disp , GetEventWeight());
1761  fhDispLam1LowE ->Fill(disp , lambda1, GetEventWeight());
1762  fhEtaLam0LowE ->Fill(eta , lambda0, GetEventWeight());
1763  fhPhiLam0LowE ->Fill(phi , lambda0, GetEventWeight());
1764  }
1765  else
1766  {
1767  fhNCellsLam0HighE ->Fill(ncells, lambda0, GetEventWeight());
1768  fhNCellsLam1HighE ->Fill(ncells, lambda1, GetEventWeight());
1769  fhNCellsDispHighE ->Fill(ncells, disp , GetEventWeight());
1770 
1771  fhLam1Lam0HighE ->Fill(lambda1, lambda0, GetEventWeight());
1772  fhLam0DispHighE ->Fill(lambda0, disp , GetEventWeight());
1773  fhDispLam1HighE ->Fill(disp , lambda1, GetEventWeight());
1774  fhEtaLam0HighE ->Fill(eta , lambda0, GetEventWeight());
1775  fhPhiLam0HighE ->Fill(phi , lambda0, GetEventWeight());
1776  }
1777  }
1778 
1779  if (IsDataMC())
1780  {
1781  AliVCaloCells* cells = 0;
1782  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
1783  else cells = GetPHOSCells();
1784 
1785  // Fill histograms to check shape of embedded clusters
1786  Float_t fraction = 0;
1787  // printf("check embedding %i\n",GetReader()->IsEmbeddedClusterSelectionOn());
1788 
1789  if(GetReader()->IsEmbeddedClusterSelectionOn())
1790  {
1791  // Only working for EMCAL
1792  // printf("embedded\n");
1793 
1794  Float_t clusterE = 0; // recalculate in case corrections applied.
1795  Float_t cellE = 0;
1796  for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
1797  {
1798  cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
1799  clusterE+=cellE;
1800  fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
1801  }
1802 
1803  //Fraction of total energy due to the embedded signal
1804  fraction/=clusterE;
1805 
1806  AliDebug(1,Form("Energy fraction of embedded signal %2.3f, Energy %2.3f",fraction, clusterE));
1807 
1808  fhEmbeddedSignalFractionEnergy->Fill(clusterE, fraction, GetEventWeight());
1809  } // embedded fraction
1810 
1811  //
1812  // Check the origin
1813  //
1814  Int_t mcIndex = -1;
1816  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
1817  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
1818  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1819  {
1820  if ( !conversion ) mcIndex = kmcssPhoton ;
1821  else mcIndex = kmcssPhotonConv;
1822 
1823  // Fill histograms to check shape of embedded clusters
1824  if(GetReader()->IsEmbeddedClusterSelectionOn())
1825  {
1826  if (fraction > 0.9)
1827  {
1828  fhEmbedPhotonELambda0FullSignal ->Fill(energy, lambda0, GetEventWeight());
1829  }
1830  else if(fraction > 0.5)
1831  {
1832  fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0, GetEventWeight());
1833  }
1834  else if(fraction > 0.1)
1835  {
1836  fhEmbedPhotonELambda0MostlyBkg ->Fill(energy, lambda0, GetEventWeight());
1837  }
1838  else
1839  {
1840  fhEmbedPhotonELambda0FullBkg ->Fill(energy, lambda0, GetEventWeight());
1841  }
1842  } // embedded
1843  } // photon no conversion
1844  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
1845  {
1846  mcIndex = kmcssElectron ;
1847  }//electron
1848  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
1849  {
1850  if ( !conversion ) mcIndex = kmcssPi0 ;
1851  else mcIndex = kmcssPi0Conv;
1852 
1853  // Fill histograms to check shape of embedded clusters
1854  if(GetReader()->IsEmbeddedClusterSelectionOn())
1855  {
1856  if (fraction > 0.9)
1857  {
1858  fhEmbedPi0ELambda0FullSignal ->Fill(energy, lambda0, GetEventWeight());
1859  }
1860  else if(fraction > 0.5)
1861  {
1862  fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0, GetEventWeight());
1863  }
1864  else if(fraction > 0.1)
1865  {
1866  fhEmbedPi0ELambda0MostlyBkg ->Fill(energy, lambda0, GetEventWeight());
1867  }
1868  else
1869  {
1870  fhEmbedPi0ELambda0FullBkg ->Fill(energy, lambda0, GetEventWeight());
1871  }
1872  } // embedded
1873 
1874  }//pi0
1875  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
1876  {
1877  if ( !conversion ) mcIndex = kmcssEta ;
1878  else mcIndex = kmcssEtaConv;
1879 
1880  }//eta
1881  else
1882  {
1883  mcIndex = kmcssOther ;
1884  }//other particles
1885 
1886  //
1887  // Fill histograms
1888  //
1889  fhMCELambda0 [mcIndex]->Fill(energy, lambda0, GetEventWeight());
1890  fhMCPtLambda0[mcIndex]->Fill(pt , lambda0, GetEventWeight());
1891  fhMCELambda1 [mcIndex]->Fill(energy, lambda1, GetEventWeight());
1892  fhMCNCellsE [mcIndex]->Fill(energy, ncells , GetEventWeight());
1893 
1895  {
1896  fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction, GetEventWeight());
1897  fhMCEDispersion [mcIndex]->Fill(energy, disp , GetEventWeight());
1898  }
1899 
1900  //
1901  // Check particle overlaps in cluster
1902  //
1903  if ( !GetReader()->IsEmbeddedClusterSelectionOn() )
1904  {
1905  // Compare the primary depositing more energy with the rest,
1906  // if no photon/electron as comon ancestor (conversions), count as other particle
1907  const UInt_t nlabels = cluster->GetNLabels();
1908  Int_t overpdg[nlabels];
1909  Int_t overlab[nlabels];
1910  Int_t noverlaps = GetMCAnalysisUtils()->GetNOverlaps(cluster->GetLabels(), nlabels,mcTag,-1,GetMC(),overpdg,overlab);
1911 
1912  //printf("N overlaps %d \n",noverlaps);
1913 
1914  if(noverlaps == 0)
1915  fhMCPtLambda0Overlaps[mcIndex][0]->Fill(pt, lambda0, GetEventWeight());
1916  else if(noverlaps == 1)
1917  fhMCPtLambda0Overlaps[mcIndex][1]->Fill(pt, lambda0, GetEventWeight());
1918  else if(noverlaps > 1)
1919  fhMCPtLambda0Overlaps[mcIndex][2]->Fill(pt, lambda0, GetEventWeight());
1920  else
1921  AliWarning(Form("n overlaps = %d!!", noverlaps));
1922  }
1923 
1925  {
1926  if (energy < 2.)
1927  {
1928  fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction, GetEventWeight());
1929  fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells, maxCellFraction, GetEventWeight());
1930  }
1931  else if(energy < 6.)
1932  {
1933  fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction, GetEventWeight());
1934  fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells, maxCellFraction, GetEventWeight());
1935  }
1936  else
1937  {
1938  fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction, GetEventWeight());
1939  fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells, maxCellFraction, GetEventWeight());
1940  }
1941 
1942  if(GetCalorimeter() == kEMCAL)
1943  {
1944  fhMCEDispEta [mcIndex]-> Fill(energy, dEta , GetEventWeight());
1945  fhMCEDispPhi [mcIndex]-> Fill(energy, dPhi , GetEventWeight());
1946  fhMCESumEtaPhi [mcIndex]-> Fill(energy, sEtaPhi , GetEventWeight());
1947  fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy, dPhi-dEta, GetEventWeight());
1948 
1949  if( dEta+dPhi > 0 ) fhMCESphericity[mcIndex]-> Fill(energy, (dPhi-dEta)/(dEta+dPhi), GetEventWeight());
1950 
1951  Int_t ebin = -1;
1952  if (energy < 2 ) ebin = 0;
1953  else if (energy < 4 ) ebin = 1;
1954  else if (energy < 6 ) ebin = 2;
1955  else if (energy < 10) ebin = 3;
1956  else if (energy < 15) ebin = 4;
1957  else if (energy < 20) ebin = 5;
1958  else ebin = 6;
1959 
1960  fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta , dPhi, GetEventWeight());
1961  fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0, dEta, GetEventWeight());
1962  fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0, dPhi, GetEventWeight());
1963  }
1964  }
1965  }// MC data
1966 }
1967 
1968 //__________________________________________________________________________
1980 //__________________________________________________________________________
1982  Int_t cut, Int_t nSMod, Bool_t matched,
1983  Int_t tag, Float_t mcbin, Float_t egen, Int_t noverlaps)
1984 {
1985  Float_t dEta = cluster->GetTrackDz();
1986  Float_t dPhi = cluster->GetTrackDx();
1987  Float_t ener = cluster->E();
1988 // if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1989 // {
1990 // dPhi = 2000., dEta = 2000.;
1991 // GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedPhiesiduals(cluster->GetID(),dEta,dPhi);
1992 // }
1993 
1994  AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1995 
1996  if ( !fhTrackMatchedDEtaPos[cut] || TMath::Abs(dPhi) > 999 || !track ) return;
1997 
1998  Bool_t positive = (track->Charge()>0);
1999  if ( positive )
2000  {
2001  fhTrackMatchedDEtaPos[cut]->Fill(ener, dEta, GetEventWeight());
2002  fhTrackMatchedDPhiPos[cut]->Fill(ener, dPhi, GetEventWeight());
2003  if(ener > 0.5) fhTrackMatchedDEtaDPhiPos[cut]->Fill(dEta,dPhi, GetEventWeight());
2004 
2005  if ( fFillTMHistoTrackPt )
2006  {
2007  fhTrackMatchedDEtaPosTrackPt[cut]->Fill(track->Pt(), dEta, GetEventWeight());
2008  fhTrackMatchedDPhiPosTrackPt[cut]->Fill(track->Pt(), dPhi, GetEventWeight());
2009  if(track->Pt() > 0.5) fhTrackMatchedDEtaDPhiPosTrackPt[cut]->Fill(dEta,dPhi, GetEventWeight());
2010  }
2011  }
2012  else
2013  {
2014  fhTrackMatchedDEtaNeg[cut]->Fill(ener, dEta, GetEventWeight());
2015  fhTrackMatchedDPhiNeg[cut]->Fill(ener, dPhi, GetEventWeight());
2016  if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg[cut]->Fill(dEta, dPhi, GetEventWeight());
2017 
2018  if ( fFillTMHistoTrackPt )
2019  {
2020  fhTrackMatchedDEtaNegTrackPt[cut]->Fill(track->Pt(), dEta, GetEventWeight());
2021  fhTrackMatchedDPhiNegTrackPt[cut]->Fill(track->Pt(), dPhi, GetEventWeight());
2022  if(track->Pt() > 0.5) fhTrackMatchedDEtaDPhiNegTrackPt[cut]->Fill(dEta, dPhi, GetEventWeight());
2023  }
2024  }
2025 
2026  // TRD
2027  if( GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
2028  nSMod >= GetFirstSMCoveredByTRD() )
2029  {
2030  fhTrackMatchedDEtaTRD[cut]->Fill(ener, dEta, GetEventWeight());
2031  fhTrackMatchedDPhiTRD[cut]->Fill(ener, dPhi, GetEventWeight());
2032  }
2033 
2034  // Check dEdx and E/p and MC origin of matched/unmatched clusters
2035  //
2036  if ( IsDataMC() )
2037  {
2038  fhTrackMatchedMCParticleBeforeTM [cut]->Fill(ener , mcbin, GetEventWeight());
2039  fhTrackMatchedMCParticleTrackPtBeforeTM[cut]->Fill(track->Pt(), mcbin, GetEventWeight());
2040  }
2041 
2042  if ( (matched && !cut) || (!matched && cut) )
2043  {
2044  Float_t dEdx = track->GetTPCsignal();
2045  Float_t eOverp = -1;
2046  if ( track->P() > 0.1 ) eOverp = ener/track->P();
2047 
2048  // if ( cut )
2049  // printf("matched %d and cut %d: E %2.2f, track p %2.2f, track pt %2.2f, dEta %2.2f, dPhi %2.2f, dEdx %2.2f, E/p %2.2f\n",
2050  // matched,cut,ener,track->P(),track->Pt(),dEta,dPhi,dEdx,eOverp);
2051 
2052  fhdEdx [cut]->Fill(ener, dEdx , GetEventWeight());
2053  fhEOverP[cut]->Fill(ener, eOverp, GetEventWeight());
2054 
2055  if ( fFillTMHistoTrackPt )
2056  {
2057  fhdEdxTrackPt [cut]->Fill(track->Pt(), dEdx , GetEventWeight());
2058  fhEOverPTrackPt[cut]->Fill(track->Pt(), eOverp, GetEventWeight());
2059  }
2060 
2061  if(GetCalorimeter()==kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
2062  nSMod >= GetFirstSMCoveredByTRD() )
2063  fhEOverPTRD[cut]->Fill(ener, eOverp, GetEventWeight());
2064 
2065  if ( IsDataMC() )
2066  {
2067  fhTrackMatchedMCParticle [cut]->Fill(ener, mcbin, GetEventWeight());
2068  fhTrackMatchedMCParticleTrackPt [cut]->Fill(track->Pt(), mcbin, GetEventWeight());
2069  fhTrackMatchedMCParticleVsEOverP[cut]->Fill(ener, mcbin, eOverp, GetEventWeight());
2070 
2071  if ( egen > 0.1 )
2072  fhTrackMatchedMCParticleVsErecEgen[cut]->Fill(ener, mcbin, ener/egen, GetEventWeight());
2073 
2075  if ( conversion )
2076  fhTrackMatchedMCParticleConverted[cut]->Fill(ener, mcbin, GetEventWeight());
2077 
2078  fhTrackMatchedMCParticleVsNOverlaps[cut]->Fill(ener, mcbin, noverlaps, GetEventWeight());
2079 
2080  if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) &&
2081  !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2082  {
2083  if ( noverlaps )
2084  {
2085  fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(ener, dEta, GetEventWeight());
2086  fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(ener, dPhi, GetEventWeight());
2087  }
2088  else
2089  {
2090  fhTrackMatchedDEtaMCOverlap[cut]->Fill(ener, dEta, GetEventWeight());
2091  fhTrackMatchedDPhiMCOverlap[cut]->Fill(ener, dPhi, GetEventWeight());
2092  }
2093 
2094  if( conversion )
2095  {
2096  fhTrackMatchedDEtaMCConversion[cut]->Fill(ener, dEta, GetEventWeight());
2097  fhTrackMatchedDPhiMCConversion[cut]->Fill(ener, dPhi, GetEventWeight());
2098  }
2099 
2100  }// Check overlaps
2101 
2102  } // MC
2103 
2104  } // track info for clusters matched and not matched
2105 
2106 }
2107 
2108 //___________________________________________
2110 //___________________________________________
2112 {
2113  TString parList ; //this will be list of parameters used for this analysis.
2114  const Int_t buffersize = 255;
2115  char onePar[buffersize] ;
2116 
2117  snprintf(onePar,buffersize,"--- AliAnaPhoton ---:") ;
2118  parList+=onePar ;
2119  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
2120  parList+=onePar ;
2121  snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster);",fMinDist) ;
2122  parList+=onePar ;
2123  snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation);",fMinDist2) ;
2124  parList+=onePar ;
2125  snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study);",fMinDist3) ;
2126  parList+=onePar ;
2127  snprintf(onePar,buffersize,"fRejectTrackMatch: %d",fRejectTrackMatch) ;
2128  parList+=onePar ;
2129 
2130  // Get parameters set in base class.
2131  parList += GetBaseParametersList() ;
2132 
2133  // Get parameters set in PID class.
2134  parList += GetCaloPID()->GetPIDParametersList() ;
2135 
2136  // Get parameters set in FiducialCut class (not available yet)
2137  //parlist += GetFidCut()->GetFidCutParametersList()
2138 
2139  return new TObjString(parList) ;
2140 }
2141 
2142 //_____________________________________________
2145 //_____________________________________________
2147 {
2148  TList * outputContainer = new TList() ;
2149  outputContainer->SetName("PhotonHistos") ;
2150 
2157 
2164 
2165  Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
2168  Int_t nPoverEbins = GetHistogramRanges()->GetHistoEOverPBins();
2169  Float_t pOverEmax = GetHistogramRanges()->GetHistoEOverPMax();
2170  Float_t pOverEmin = GetHistogramRanges()->GetHistoEOverPMin();
2171 
2172  Int_t nratbins = GetHistogramRanges()->GetHistoRatioBins();
2175  Int_t ndifbins = GetHistogramRanges()->GetHistoEDiffBins();
2178 
2182 
2186 
2187  // Init the number of modules, set in the class AliCalorimeterUtils
2188  //
2189  InitCaloParameters(); // See AliCaloTrackCorrBaseClass
2190 
2191  Int_t totalSM = fLastModule-fFirstModule+1;
2192 
2193  //printf("N SM %d, first SM %d, last SM %d, total %d\n",fNModules,fFirstModule,fLastModule, totalSM);
2194 
2195  // Cell column-row histograms, see base class for data members setting
2196  //fNMaxColsFull+2,-1.5,fNMaxColsFull+0.5, fNMaxRowsFull+2,-1.5,fNMaxRowsFull+0.5
2197  Int_t ncolcell = fNMaxColsFull+2;
2198  Float_t colcellmin = -1.5;
2199  Float_t colcellmax = fNMaxColsFull+0.5;
2200 
2202  Float_t rowcellmin = fNMaxRowsFullMin-1.5;
2203  Float_t rowcellmax = fNMaxRowsFullMax+0.5;
2204 
2205  Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins for SS studies
2206 
2207  TString cut[] = {"Open","Reader","E","Time","NCells","NLM","Fidutial","Matching","Bad","PID"};
2208  for (Int_t i = 0; i < 10 ; i++)
2209  {
2210  fhClusterCutsE[i] = new TH1F(Form("hE_Cut_%d_%s", i, cut[i].Data()),
2211  Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
2212  nptbins,ptmin,ptmax);
2213  fhClusterCutsE[i]->SetYTitle("d#it{N}/d#it{E} ");
2214  fhClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
2215  outputContainer->Add(fhClusterCutsE[i]) ;
2216 
2217  fhClusterCutsPt[i] = new TH1F(Form("hPt_Cut_%d_%s", i, cut[i].Data()),
2218  Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
2219  nptbins,ptmin,ptmax);
2220  fhClusterCutsPt[i]->SetYTitle("d#it{N}/d#it{E} ");
2221  fhClusterCutsPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2222  outputContainer->Add(fhClusterCutsPt[i]) ;
2223  }
2224 
2225  fhEClusterSM = new TH2F("hEClusterSM","Raw clusters E and super-module number",
2226  nptbins,ptmin,ptmax,
2227  totalSM,fFirstModule-0.5,fLastModule+0.5);
2228  fhEClusterSM->SetYTitle("SuperModule ");
2229  fhEClusterSM->SetXTitle("#it{E} (GeV)");
2230  outputContainer->Add(fhEClusterSM) ;
2231 
2232  fhPtClusterSM = new TH2F("hPtClusterSM","Raw clusters #it{p}_{T} and super-module number",
2233  nptbins,ptmin,ptmax,
2234  totalSM,fFirstModule-0.5,fLastModule+0.5);
2235  fhPtClusterSM->SetYTitle("SuperModule ");
2236  fhPtClusterSM->SetXTitle("#it{E} (GeV)");
2237  outputContainer->Add(fhPtClusterSM) ;
2238 
2239  fhEPhotonSM = new TH2F("hEPhotonSM","Selected clusters E and super-module number",
2240  nptbins,ptmin,ptmax,
2241  totalSM,fFirstModule-0.5,fLastModule+0.5);
2242  fhEPhotonSM->SetYTitle("SuperModule ");
2243  fhEPhotonSM->SetXTitle("#it{E} (GeV)");
2244  outputContainer->Add(fhEPhotonSM) ;
2245 
2246  fhPtPhotonSM = new TH2F("hPtPhotonSM","Selected clusters #it{p}_{T} and super-module number",
2247  nptbins,ptmin,ptmax,
2248  totalSM,fFirstModule-0.5,fLastModule+0.5);
2249  fhPtPhotonSM->SetYTitle("SuperModule ");
2250  fhPtPhotonSM->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2251  outputContainer->Add(fhPtPhotonSM) ;
2252 
2254  {
2255 // printf("Bins Mult: n %d, min %d, max %d; Sum: n %d, min %f, max %f\n",
2256 // nmultbin,multmin,multmax,nsumbin,summin, summax);
2257  for(Int_t icut = 0; icut < GetReader()->GetTrackMultiplicityNPtCut(); icut++)
2258  {
2259  fhPtPhotonNTracks[icut] = new TH2F
2260  (Form("hPtPhotonNTracks_PtCut%d",icut),
2261  Form("Number of tracks per event with |#eta|<%2.2f and #it{p}_{T} > %2.2f GeV/#it{c}",
2262  GetReader()->GetTrackMultiplicityEtaCut(),GetReader()->GetTrackMultiplicityPtCut(icut)),
2263  nptbins,ptmin,ptmax, nmultbin,multmin,multmax);
2264  fhPtPhotonNTracks[icut]->SetYTitle("# of tracks");
2265  fhPtPhotonNTracks[icut]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2266  outputContainer->Add(fhPtPhotonNTracks[icut]);
2267 
2268  fhPtPhotonSumPtTracks[icut] = new TH2F
2269  (Form("hPtPhotonSumPtTracks_PtCut%d",icut),
2270  Form("#Sigma #it{p}_{T} of tracks per event with |#eta|<%2.2f and #it{p}_{T} > %2.2f GeV/#it{c}",
2271  GetReader()->GetTrackMultiplicityEtaCut(),GetReader()->GetTrackMultiplicityPtCut(icut)),
2272  nptbins,ptmin,ptmax, nsumbin,summin,summax);
2273  fhPtPhotonSumPtTracks[icut]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2274  fhPtPhotonSumPtTracks[icut]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2275  outputContainer->Add(fhPtPhotonSumPtTracks[icut]);
2276  }
2277  }
2278 
2279  fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
2280  fhNCellsE->SetXTitle("#it{E} (GeV)");
2281  fhNCellsE->SetYTitle("# of cells in cluster");
2282  outputContainer->Add(fhNCellsE);
2283 
2284  fhCellsE = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
2285  fhCellsE->SetXTitle("#it{E}_{cluster} (GeV)");
2286  fhCellsE->SetYTitle("#it{E}_{cell} (GeV)");
2287  outputContainer->Add(fhCellsE);
2288 
2289  fhTimePt = new TH2F ("hTimePt","time of cluster vs pT of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2290  fhTimePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2291  fhTimePt->SetYTitle("#it{time} (ns)");
2292  outputContainer->Add(fhTimePt);
2293 
2295  {
2296  fhMaxCellDiffClusterE = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
2297  nptbins,ptmin,ptmax, 500,0,1.);
2298  fhMaxCellDiffClusterE->SetXTitle("#it{E}_{cluster} (GeV) ");
2299  fhMaxCellDiffClusterE->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
2300  outputContainer->Add(fhMaxCellDiffClusterE);
2301  }
2302 
2303  fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
2304  fhEPhoton->SetYTitle("#it{counts}");
2305  fhEPhoton->SetXTitle("#it{E}_{#gamma}(GeV)");
2306  outputContainer->Add(fhEPhoton) ;
2307 
2308  fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs #it{p}_{T}",nptbins,ptmin,ptmax);
2309  fhPtPhoton->SetYTitle("#it{counts}");
2310  fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/#it{c})");
2311  outputContainer->Add(fhPtPhoton) ;
2312 
2314  {
2315  fhPtCentralityPhoton = new TH2F("hPtCentralityPhoton","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
2316  fhPtCentralityPhoton->SetYTitle("Centrality");
2317  fhPtCentralityPhoton->SetXTitle("#it{p}_{T}(GeV/#it{c})");
2318  outputContainer->Add(fhPtCentralityPhoton) ;
2319 
2320  fhPtEventPlanePhoton = new TH2F("hPtEventPlanePhoton","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
2321  fhPtEventPlanePhoton->SetYTitle("Event plane angle (rad)");
2322  fhPtEventPlanePhoton->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2323  outputContainer->Add(fhPtEventPlanePhoton) ;
2324  }
2325 
2326  fhEtaPhi = new TH2F
2327  ("hEtaPhi","cluster,#it{E} > 0.5 GeV, #eta vs #varphi",netabins,etamin,etamax,nphibins,phimin,phimax);
2328  fhEtaPhi->SetYTitle("#varphi (rad)");
2329  fhEtaPhi->SetXTitle("#eta");
2330  outputContainer->Add(fhEtaPhi) ;
2331 
2332  fhPhiPhoton = new TH2F
2333  ("hPhiPhoton","#varphi_{#gamma} vs #it{p}_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2334  fhPhiPhoton->SetYTitle("#varphi (rad)");
2335  fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
2336  outputContainer->Add(fhPhiPhoton) ;
2337 
2338  fhEtaPhoton = new TH2F
2339  ("hEtaPhoton","#eta_{#gamma} vs #it{p}_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2340  fhEtaPhoton->SetYTitle("#eta");
2341  fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
2342  outputContainer->Add(fhEtaPhoton) ;
2343 
2344  fhEtaPhiPhoton = new TH2F
2345  ("hEtaPhiPhoton","#eta vs #varphi",netabins,etamin,etamax,nphibins,phimin,phimax);
2346  fhEtaPhiPhoton->SetYTitle("#varphi (rad)");
2347  fhEtaPhiPhoton->SetXTitle("#eta");
2348  outputContainer->Add(fhEtaPhiPhoton) ;
2349  if(GetMinPt() < 0.5)
2350  {
2351  fhEtaPhi05Photon = new TH2F
2352  ("hEtaPhi05Photon","#eta vs #varphi, E < 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
2353  fhEtaPhi05Photon->SetYTitle("#varphi (rad)");
2354  fhEtaPhi05Photon->SetXTitle("#eta");
2355  outputContainer->Add(fhEtaPhi05Photon) ;
2356  }
2357 
2358  fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
2359  nptbins,ptmin,ptmax,10,0,10);
2360  fhNLocMax ->SetYTitle("N maxima");
2361  fhNLocMax ->SetXTitle("#it{E} (GeV)");
2362  outputContainer->Add(fhNLocMax) ;
2363 
2364  TString mcPartLabels[] = {"#gamma","#pi^{0}, #eta (#gamma#gamma)","e^{#pm}","#mu^{#pm}",
2365  "#pi^{#pm}","k^{#pm}","p","#bar{p}","n","#bar{n}","Other"};
2366  TString cutCase[] = {"Raw","BeforeTM","AfterTM","Final"};
2367  if ( IsDataMC() )
2368  {
2369  for(Int_t i = 0; i < 4; i++)
2370  {
2371  fhMCParticle[i] = new TH2F
2372  (Form("hMCParticle%s",cutCase[i].Data()),
2373  Form("Origin of particle vs energy %s",cutCase[i].Data()),
2374  nptbins,ptmin,ptmax,11,0,11);
2375  fhMCParticle[i]->SetXTitle("#it{E} (GeV)");
2376  //fhMCParticle[i]->SetYTitle("Particle type");
2377  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
2378  fhMCParticle[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
2379 
2380  outputContainer->Add(fhMCParticle[i]);
2381 
2382  fhMCParticleConverted[i] = new TH2F
2383  (Form("hMCParticle%s_Converted",cutCase[i].Data()),
2384  Form("Origin of particle vs energy %s, suffer conversion",cutCase[i].Data()),
2385  nptbins,ptmin,ptmax,11,0,11);
2386  fhMCParticleConverted[i]->SetXTitle("#it{E} (GeV)");
2387  //fhMCParticle[i]->SetYTitle("Particle type");
2388  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
2389  fhMCParticleConverted[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
2390 
2391  outputContainer->Add(fhMCParticleConverted[i]);
2392 
2393  fhMCParticleVsNOverlaps[i] = new TH3F
2394  (Form("hMCParticle%s_NOverlaps",cutCase[i].Data()),
2395  Form("Origin of particle vs energy %s vs number of particle overlaps",cutCase[i].Data()),
2396  nptbins,ptmin,ptmax,11,0,11,10,0,10);
2397  fhMCParticleVsNOverlaps[i]->SetXTitle("#it{E} (GeV)");
2398  fhMCParticleVsNOverlaps[i]->SetZTitle("#it{n}_{overlaps}");
2399  //fhMCParticle[i]->SetYTitle("Particle type");
2400  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
2401  fhMCParticleVsNOverlaps[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
2402 
2403  outputContainer->Add(fhMCParticleVsNOverlaps[i]);
2404 
2405  fhMCParticleVsErecEgen[i] = new TH3F
2406  (Form("hMCParticle%s_ErecEgenRat",cutCase[i].Data()),
2407  Form("Origin of particle vs energy %s vs #it{E}_{rec}/#it{E}_{gen}",cutCase[i].Data()),
2408  nptbins,ptmin,ptmax,11,0,11,100,0,2);
2409  fhMCParticleVsErecEgen[i]->SetXTitle("#it{E}_{rec} (GeV)");
2410  fhMCParticleVsErecEgen[i]->SetZTitle("#it{E}_{rec}/#it{E}_{gen}");
2411  //fhMCParticle[i]->SetYTitle("Particle type");
2412  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
2413  fhMCParticleVsErecEgen[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
2414 
2415  outputContainer->Add(fhMCParticleVsErecEgen[i]);
2416  }
2417  }
2418 
2419 
2420  //Shower shape
2421  if(fFillSSHistograms)
2422  {
2423  fhLam0E = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2424  fhLam0E->SetYTitle("#lambda_{0}^{2}");
2425  fhLam0E->SetXTitle("#it{E} (GeV)");
2426  outputContainer->Add(fhLam0E);
2427 
2428  fhLam0Pt = new TH2F ("hLam0Pt","#lambda_{0}^{2} vs #it{p}_{T}", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2429  fhLam0Pt->SetYTitle("#lambda_{0}^{2}");
2430  fhLam0Pt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2431  outputContainer->Add(fhLam0Pt);
2432 
2433  fhLam1E = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2434  fhLam1E->SetYTitle("#lambda_{1}^{2}");
2435  fhLam1E->SetXTitle("#it{E} (GeV)");
2436  outputContainer->Add(fhLam1E);
2437 
2438  fhLam1Pt = new TH2F ("hLam1Pt","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2439  fhLam1Pt->SetYTitle("#lambda_{1}^{2}");
2440  fhLam1Pt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2441  outputContainer->Add(fhLam1Pt);
2442 
2444  {
2445  for(Int_t ism = 0; ism < fNModules; ism++)
2446  {
2447  if(ism < fFirstModule || ism > fLastModule) continue;
2448  fhLam0PerSM[ism] = new TH2F
2449  (Form("hLam0_SM%d",ism),
2450  Form("#it{p}_{T} vs #lambda^{2}_{0} in SM %d",ism),
2451  nptbins,ptmin,ptmax,40,0,0.4);
2452  fhLam0PerSM[ism]->SetYTitle("#lambda^{2}_{0}");
2453  fhLam0PerSM[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2454  outputContainer->Add(fhLam0PerSM[ism]) ;
2455 
2456  fhLam1PerSM[ism] = new TH2F
2457  (Form("hLam1_SM%d",ism),
2458  Form("#it{p}_{T} vs #lambda^{2}_{1} in SM %d",ism),
2459  nptbins,ptmin,ptmax,40,0,0.4);
2460  fhLam1PerSM[ism]->SetYTitle("#lambda^{2}_{1}");
2461  fhLam1PerSM[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2462  outputContainer->Add(fhLam1PerSM[ism]) ;
2463  }
2464  }
2465 
2467  {
2468  fhDispE = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2469  fhDispE->SetYTitle("D^{2}");
2470  fhDispE->SetXTitle("#it{E} (GeV) ");
2471  outputContainer->Add(fhDispE);
2472 
2473  fhDispPt = new TH2F ("hDispPt"," dispersion^{2} vs #it{p}_{T}", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2474  fhDispPt->SetYTitle("D^{2}");
2475  fhDispPt->SetXTitle("#it{p}_{T} (GeV/#it{c}) ");
2476  outputContainer->Add(fhDispPt);
2477  }
2478 
2480  {
2481  fhLam0PtNLM1 = new TH2F ("hLam0PtNLM1","#lambda_{0}^{2} vs #it{p}_{T}, #it{n}_{LM}=1", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2482  fhLam0PtNLM1->SetYTitle("#lambda_{0}^{2}");
2483  fhLam0PtNLM1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2484  outputContainer->Add(fhLam0PtNLM1);
2485 
2486  fhLam0PtNLM2 = new TH2F ("hLam0PtNLM2","#lambda_{0}^{2} vs #it{p}_{T}, #it{n}_{LM}=2", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2487  fhLam0PtNLM2->SetYTitle("#lambda_{0}^{2}");
2488  fhLam0PtNLM2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2489  outputContainer->Add(fhLam0PtNLM2);
2490 
2491  fhLam1PtNLM1 = new TH2F ("hLam1PtNLM1","#lambda_{1}^{2} vs #it{p}_{T}, #it{n}_{LM}=1", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2492  fhLam1PtNLM1->SetYTitle("#lambda_{1}^{2}");
2493  fhLam1PtNLM1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2494  outputContainer->Add(fhLam1PtNLM1);
2495 
2496  fhLam1PtNLM2 = new TH2F ("hLam1PtNLM2","#lambda_{1}^{2} vs #it{p}_{T}, #it{n}_{LM}=2", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2497  fhLam1PtNLM2->SetYTitle("#lambda_{1}^{2}");
2498  fhLam1PtNLM2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2499  outputContainer->Add(fhLam1PtNLM2);
2500  }
2501 
2502  if(!fRejectTrackMatch)
2503  {
2504  fhLam0ETM = new TH2F ("hLam0ETM","#lambda_{0}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2505  fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
2506  fhLam0ETM->SetXTitle("#it{E} (GeV)");
2507  outputContainer->Add(fhLam0ETM);
2508 
2509  fhLam0PtTM = new TH2F ("hLam0PtTM","#lambda_{0}^{2} vs #it{p}_{T}, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2510  fhLam0PtTM->SetYTitle("#lambda_{0}^{2}");
2511  fhLam0PtTM->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2512  outputContainer->Add(fhLam0PtTM);
2513 
2514  fhLam1ETM = new TH2F ("hLam1ETM","#lambda_{1}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2515  fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
2516  fhLam1ETM->SetXTitle("#it{E} (GeV)");
2517  outputContainer->Add(fhLam1ETM);
2518 
2520  {
2521  fhDispETM = new TH2F ("hDispETM"," dispersion^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2522  fhDispETM->SetYTitle("D^{2}");
2523  fhDispETM->SetXTitle("#it{E} (GeV) ");
2524  outputContainer->Add(fhDispETM);
2525  }
2526  }
2527 
2528  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
2529  {
2530  fhLam0ETRD = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2531  fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
2532  fhLam0ETRD->SetXTitle("#it{E} (GeV)");
2533  outputContainer->Add(fhLam0ETRD);
2534 
2535  fhLam0PtTRD = new TH2F ("hLam0PtTRD","#lambda_{0}^{2} vs #it{p}_{T}, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2536  fhLam0PtTRD->SetYTitle("#lambda_{0}^{2}");
2537  fhLam0PtTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2538  outputContainer->Add(fhLam0PtTRD);
2539 
2540  fhLam1ETRD = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2541  fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
2542  fhLam1ETRD->SetXTitle("#it{E} (GeV)");
2543  outputContainer->Add(fhLam1ETRD);
2544 
2546  {
2547  fhDispETRD = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2548  fhDispETRD->SetYTitle("Dispersion^{2}");
2549  fhDispETRD->SetXTitle("#it{E} (GeV) ");
2550  outputContainer->Add(fhDispETRD);
2551  }
2552 
2554  {
2555  fhLam0ETMTRD = new TH2F ("hLam0ETMTRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2556  fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
2557  fhLam0ETMTRD->SetXTitle("#it{E} (GeV)");
2558  outputContainer->Add(fhLam0ETMTRD);
2559 
2560  fhLam0PtTMTRD = new TH2F ("hLam0PtTMTRD","#lambda_{0}^{2} vs #it{p}_{T}, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2561  fhLam0PtTMTRD->SetYTitle("#lambda_{0}^{2}");
2562  fhLam0PtTMTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2563  outputContainer->Add(fhLam0PtTMTRD);
2564 
2565  fhLam1ETMTRD = new TH2F ("hLam1ETMTRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2566  fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
2567  fhLam1ETMTRD->SetXTitle("#it{E} (GeV)");
2568  outputContainer->Add(fhLam1ETMTRD);
2569 
2571  {
2572  fhDispETMTRD = new TH2F ("hDispETMTRD"," dispersion^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #varphi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2573  fhDispETMTRD->SetYTitle("Dispersion^{2}");
2574  fhDispETMTRD->SetXTitle("#it{E} (GeV) ");
2575  outputContainer->Add(fhDispETMTRD);
2576  }
2577  }
2578  }
2579 
2581  {
2582  fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2583  fhNCellsLam0LowE->SetXTitle("N_{cells}");
2584  fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
2585  outputContainer->Add(fhNCellsLam0LowE);
2586 
2587  fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2588  fhNCellsLam0HighE->SetXTitle("N_{cells}");
2589  fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
2590  outputContainer->Add(fhNCellsLam0HighE);
2591 
2592  fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2593  fhNCellsLam1LowE->SetXTitle("N_{cells}");
2594  fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
2595  outputContainer->Add(fhNCellsLam1LowE);
2596 
2597  fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2598  fhNCellsLam1HighE->SetXTitle("N_{cells}");
2599  fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
2600  outputContainer->Add(fhNCellsLam1HighE);
2601 
2602  fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2603  fhNCellsDispLowE->SetXTitle("N_{cells}");
2604  fhNCellsDispLowE->SetYTitle("D^{2}");
2605  outputContainer->Add(fhNCellsDispLowE);
2606 
2607  fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2608  fhNCellsDispHighE->SetXTitle("N_{cells}");
2609  fhNCellsDispHighE->SetYTitle("D^{2}");
2610  outputContainer->Add(fhNCellsDispHighE);
2611 
2612  fhEtaLam0LowE = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
2613  fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
2614  fhEtaLam0LowE->SetXTitle("#eta");
2615  outputContainer->Add(fhEtaLam0LowE);
2616 
2617  fhPhiLam0LowE = new TH2F ("hPhiLam0LowE","#varphi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
2618  fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
2619  fhPhiLam0LowE->SetXTitle("#varphi (rad)");
2620  outputContainer->Add(fhPhiLam0LowE);
2621 
2622  fhEtaLam0HighE = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, #it{E} > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
2623  fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
2624  fhEtaLam0HighE->SetXTitle("#eta");
2625  outputContainer->Add(fhEtaLam0HighE);
2626 
2627  fhPhiLam0HighE = new TH2F ("hPhiLam0HighE","#varphi vs #lambda_{0}^{2}, #it{E} > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
2628  fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
2629  fhPhiLam0HighE->SetXTitle("#varphi (rad)");
2630  outputContainer->Add(fhPhiLam0HighE);
2631 
2632  fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2633  fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
2634  fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
2635  outputContainer->Add(fhLam1Lam0LowE);
2636 
2637  fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2638  fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
2639  fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
2640  outputContainer->Add(fhLam1Lam0HighE);
2641 
2642  fhLam0DispLowE = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2643  fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
2644  fhLam0DispLowE->SetYTitle("D^{2}");
2645  outputContainer->Add(fhLam0DispLowE);
2646 
2647  fhLam0DispHighE = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2648  fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
2649  fhLam0DispHighE->SetYTitle("D^{2}");
2650  outputContainer->Add(fhLam0DispHighE);
2651 
2652  fhDispLam1LowE = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2653  fhDispLam1LowE->SetXTitle("D^{2}");
2654  fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
2655  outputContainer->Add(fhDispLam1LowE);
2656 
2657  fhDispLam1HighE = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2658  fhDispLam1HighE->SetXTitle("D^{2}");
2659  fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
2660  outputContainer->Add(fhDispLam1HighE);
2661 
2662  if(GetCalorimeter() == kEMCAL)
2663  {
2664  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);
2665  fhDispEtaE->SetXTitle("#it{E} (GeV)");
2666  fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
2667  outputContainer->Add(fhDispEtaE);
2668 
2669  fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i} - <#varphi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2670  fhDispPhiE->SetXTitle("#it{E} (GeV)");
2671  fhDispPhiE->SetYTitle("#sigma^{2}_{#varphi #varphi}");
2672  outputContainer->Add(fhDispPhiE);
2673 
2674  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);
2675  fhSumEtaE->SetXTitle("#it{E} (GeV)");
2676  fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
2677  outputContainer->Add(fhSumEtaE);
2678 
2679  fhSumPhiE = new TH2F ("hSumPhiE","#delta^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i})^{2}/ #Sigma w_{i} - <#varphi>^{2} vs E",
2680  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2681  fhSumPhiE->SetXTitle("#it{E} (GeV)");
2682  fhSumPhiE->SetYTitle("#delta^{2}_{#varphi #varphi}");
2683  outputContainer->Add(fhSumPhiE);
2684 
2685  fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #varphi} = #Sigma w_{i}(#varphi_{i} #eta_{i} ) / #Sigma w_{i} - <#varphi><#eta> vs E",
2686  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
2687  fhSumEtaPhiE->SetXTitle("#it{E} (GeV)");
2688  fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #varphi}");
2689  outputContainer->Add(fhSumEtaPhiE);
2690 
2691  fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta} vs E",
2692  nptbins,ptmin,ptmax,200, -10,10);
2693  fhDispEtaPhiDiffE->SetXTitle("#it{E} (GeV)");
2694  fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#varphi #varphi}-#sigma^{2}_{#eta #eta}");
2695  outputContainer->Add(fhDispEtaPhiDiffE);
2696 
2697  fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi}) vs E",
2698  nptbins,ptmin,ptmax, 200, -1,1);
2699  fhSphericityE->SetXTitle("#it{E} (GeV)");
2700  fhSphericityE->SetYTitle("s = (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi})");
2701  outputContainer->Add(fhSphericityE);
2702 
2703  fhDispSumEtaDiffE = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
2704  fhDispSumEtaDiffE->SetXTitle("#it{E} (GeV)");
2705  fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average");
2706  outputContainer->Add(fhDispSumEtaDiffE);
2707 
2708  fhDispSumPhiDiffE = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#varphi #varphi} - #delta^{2}_{#varphi #varphi} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
2709  fhDispSumPhiDiffE->SetXTitle("#it{E} (GeV)");
2710  fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#varphi #varphi} - #delta^{2}_{#varphi #varphi} / average");
2711  outputContainer->Add(fhDispSumPhiDiffE);
2712 
2713  for(Int_t i = 0; i < 7; i++)
2714  {
2715  fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#varphi #varphi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
2716  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2717  fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2718  fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
2719  outputContainer->Add(fhDispEtaDispPhi[i]);
2720 
2721  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]),
2722  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2723  fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
2724  fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
2725  outputContainer->Add(fhLambda0DispEta[i]);
2726 
2727  fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#varphi #varphi} for %d < E < %d GeV",bin[i],bin[i+1]),
2728  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2729  fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
2730  fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
2731  outputContainer->Add(fhLambda0DispPhi[i]);
2732  }
2733  }
2734  }
2735  } // Shower shape
2736 
2737  // Track Matching
2738 
2739  if ( fFillTMHisto )
2740  {
2741  TString cutTM [] = {"NoCut",""};
2742 
2743  for(Int_t i = 0; i < 2; i++)
2744  {
2745 // fhTrackMatchedDEta[i] = new TH2F
2746 // (Form("hTrackMatchedDEta%s",cutTM[i].Data()),
2747 // Form("d#eta of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2748 // nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2749 // fhTrackMatchedDEta[i]->SetYTitle("#Delta #eta");
2750 // fhTrackMatchedDEta[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2751 //
2752 // fhTrackMatchedDPhi[i] = new TH2F
2753 // (Form("hTrackMatchedDPhi%s",cutTM[i].Data()),
2754 // Form("#Delta #varphi of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2755 // nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2756 // fhTrackMatchedDPhi[i]->SetYTitle("#Delta #varphi (rad)");
2757 // fhTrackMatchedDPhi[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2758 //
2759 // fhTrackMatchedDEtaDPhi[i] = new TH2F
2760 // (Form("hTrackMatchedDEtaDPhi%s",cutTM[i].Data()),
2761 // Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{E}_{cluster}, %s",cutTM[i].Data()),
2762 // nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2763 // fhTrackMatchedDEtaDPhi[i]->SetYTitle("#Delta #varphi (rad)");
2764 // fhTrackMatchedDEtaDPhi[i]->SetXTitle("#Delta #eta");
2765 
2766  fhTrackMatchedDEtaPos[i] = new TH2F
2767  (Form("hTrackMatchedDEtaPos%s",cutTM[i].Data()),
2768  Form("#Delta #eta of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2769  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2770  fhTrackMatchedDEtaPos[i]->SetYTitle("#Delta #eta");
2771  fhTrackMatchedDEtaPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2772 
2773  fhTrackMatchedDPhiPos[i] = new TH2F
2774  (Form("hTrackMatchedDPhiPos%s",cutTM[i].Data()),
2775  Form("#Delta #varphi of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2776  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2777  fhTrackMatchedDPhiPos[i]->SetYTitle("#Delta #varphi (rad)");
2778  fhTrackMatchedDPhiPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2779 
2781  (Form("hTrackMatchedDEtaDPhiPos%s",cutTM[i].Data()),
2782  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{E}_{cluster}, %s",cutTM[i].Data()),
2783  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2784  fhTrackMatchedDEtaDPhiPos[i]->SetYTitle("#Delta #varphi (rad)");
2785  fhTrackMatchedDEtaDPhiPos[i]->SetXTitle("#Delta #eta");
2786 
2787  fhTrackMatchedDEtaNeg[i] = new TH2F
2788  (Form("hTrackMatchedDEtaNeg%s",cutTM[i].Data()),
2789  Form("#Delta #eta of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2790  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2791  fhTrackMatchedDEtaNeg[i]->SetYTitle("#Delta #eta");
2792  fhTrackMatchedDEtaNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2793 
2794  fhTrackMatchedDPhiNeg[i] = new TH2F
2795  (Form("hTrackMatchedDPhiNeg%s",cutTM[i].Data()),
2796  Form("#Delta #varphi of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2797  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2798  fhTrackMatchedDPhiNeg[i]->SetYTitle("#Delta #varphi (rad)");
2799  fhTrackMatchedDPhiNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2800 
2802  (Form("hTrackMatchedDEtaDPhiNeg%s",cutTM[i].Data()),
2803  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{E}_{cluster} > 0.5 GeV, %s",cutTM[i].Data()),
2804  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2805  fhTrackMatchedDEtaDPhiNeg[i]->SetYTitle("#Delta #varphi (rad)");
2806  fhTrackMatchedDEtaDPhiNeg[i]->SetXTitle("#Delta #eta");
2807 
2808  fhdEdx[i] = new TH2F (Form("hdEdx%s",cutTM[i].Data()),Form("matched track <d#it{E}/d#it{x}> vs cluster #it{E}, %s",cutTM[i].Data()),
2809  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2810  fhdEdx[i]->SetXTitle("#it{E}^{cluster} (GeV)");
2811  fhdEdx[i]->SetYTitle("<d#it{E}/d#it{x}>");
2812 
2813  fhEOverP[i] = new TH2F (Form("hEOverP%s",cutTM[i].Data()),Form("matched track #it{E}/#it{p} vs cluster #it{E}, %s",cutTM[i].Data()),
2814  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2815  fhEOverP[i]->SetXTitle("#it{E}^{cluster} (GeV)");
2816  fhEOverP[i]->SetYTitle("#it{E}/#it{p}");
2817 
2818 // outputContainer->Add(fhTrackMatchedDEta[i]) ;
2819 // outputContainer->Add(fhTrackMatchedDPhi[i]) ;
2820 // outputContainer->Add(fhTrackMatchedDEtaDPhi[i]) ;
2821  outputContainer->Add(fhTrackMatchedDEtaPos[i]) ;
2822  outputContainer->Add(fhTrackMatchedDPhiPos[i]) ;
2823  outputContainer->Add(fhTrackMatchedDEtaDPhiPos[i]) ;
2824  outputContainer->Add(fhTrackMatchedDEtaNeg[i]) ;
2825  outputContainer->Add(fhTrackMatchedDPhiNeg[i]) ;
2826  outputContainer->Add(fhTrackMatchedDEtaDPhiNeg[i]) ;
2827  outputContainer->Add(fhdEdx[i]);
2828  outputContainer->Add(fhEOverP[i]);
2829 
2831  {
2832 // fhTrackMatchedDEtaTrackPt[i] = new TH2F
2833 // (Form("hTrackMatchedDEtaTrackPt%s",cutTM[i].Data()),
2834 // Form("#Delta #eta of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2835 // nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2836 // fhTrackMatchedDEtaTrackPt[i]->SetYTitle("#Delta #eta");
2837 // fhTrackMatchedDEtaTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2838 //
2839 // fhTrackMatchedDPhiTrackPt[i] = new TH2F
2840 // (Form("hTrackMatchedDPhiTrackPt%s",cutTM[i].Data()),
2841 // Form("#Delta #varphi of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2842 // nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2843 // fhTrackMatchedDPhiTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2844 // fhTrackMatchedDPhiTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2845 //
2846 // fhTrackMatchedDEtaDPhiTrackPt[i] = new TH2F
2847 // (Form("hTrackMatchedDEtaDPhiTrackPt%s",cutTM[i].Data()),
2848 // Form("#Delta #eta vs #Delta #varphi of cluster-track #it{p}_{T}^{track} > 0.5 GeV/#it{c}, %s",cutTM[i].Data()),
2849 // nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2850 // fhTrackMatchedDEtaDPhiTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2851 // fhTrackMatchedDEtaDPhiTrackPt[i]->SetXTitle("#Delta #eta");
2852 
2854  (Form("hTrackMatchedDEtaPosTrackPt%s",cutTM[i].Data()),
2855  Form("#Delta #eta of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2856  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2857  fhTrackMatchedDEtaPosTrackPt[i]->SetYTitle("#Delta #eta");
2858  fhTrackMatchedDEtaPosTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2859 
2861  (Form("hTrackMatchedDPhiPosTrackPt%s",cutTM[i].Data()),
2862  Form("#Delta #varphi of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2863  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2864  fhTrackMatchedDPhiPosTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2865  fhTrackMatchedDPhiPosTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2866 
2868  (Form("hTrackMatchedDEtaDPhiPosTrackPt%s",cutTM[i].Data()),
2869  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{p}_{T}^{track} > 0.5 GeV/#it{c}, %s",cutTM[i].Data()),
2870  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2871  fhTrackMatchedDEtaDPhiPosTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2872  fhTrackMatchedDEtaDPhiPosTrackPt[i]->SetXTitle("#Delta #eta");
2873 
2875  (Form("hTrackMatchedDEtaNegTrackPt%s",cutTM[i].Data()),
2876  Form("#Delta #eta of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2877  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2878  fhTrackMatchedDEtaNegTrackPt[i]->SetYTitle("#Delta #eta");
2879  fhTrackMatchedDEtaNegTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2880 
2882  (Form("hTrackMatchedDPhiNegTrackPt%s",cutTM[i].Data()),
2883  Form("#Delta #varphi of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2884  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2885  fhTrackMatchedDPhiNegTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2886  fhTrackMatchedDPhiNegTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2887 
2889  (Form("hTrackMatchedDEtaDPhiNegTrackPt%s",cutTM[i].Data()),
2890  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{p}_{T}^{track} > 0.5 GeV/#it{c}, %s",cutTM[i].Data()),
2891  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2892  fhTrackMatchedDEtaDPhiNegTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2893  fhTrackMatchedDEtaDPhiNegTrackPt[i]->SetXTitle("#Delta #eta");
2894 
2895  fhdEdxTrackPt[i] = new TH2F (Form("hdEdxTrackPt%s",cutTM[i].Data()),Form("matched track <d#it{E}/d#it{x}> vs track #it{p}_{T}, %s",cutTM[i].Data()),
2896  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2897  fhdEdxTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2898  fhdEdxTrackPt[i]->SetYTitle("<d#it{E}/d#it{x}>");
2899 
2900  fhEOverPTrackPt[i] = new TH2F (Form("hEOverPTrackPt%s",cutTM[i].Data()),Form("matched track #it{E}/#it{p} vs track #it{p}_{T}, %s",cutTM[i].Data()),
2901  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2902  fhEOverPTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2903  fhEOverPTrackPt[i]->SetYTitle("#it{E}/#it{p}");
2904 
2905 // outputContainer->Add(fhTrackMatchedDEtaTrackPt[i]) ;
2906 // outputContainer->Add(fhTrackMatchedDPhiTrackPt[i]) ;
2907 // outputContainer->Add(fhTrackMatchedDEtaDPhiTrackPt[i]) ;
2908  outputContainer->Add(fhTrackMatchedDEtaPosTrackPt[i]) ;
2909  outputContainer->Add(fhTrackMatchedDPhiPosTrackPt[i]) ;
2910  outputContainer->Add(fhTrackMatchedDEtaDPhiPosTrackPt[i]) ;
2911  outputContainer->Add(fhTrackMatchedDEtaNegTrackPt[i]) ;
2912  outputContainer->Add(fhTrackMatchedDPhiNegTrackPt[i]) ;
2913  outputContainer->Add(fhTrackMatchedDEtaDPhiNegTrackPt[i]) ;
2914  outputContainer->Add(fhdEdxTrackPt[i]);
2915  outputContainer->Add(fhEOverPTrackPt[i]);
2916  }
2917 
2919  {
2920  fhTrackMatchedDEtaTRD[i] = new TH2F
2921  (Form("hTrackMatchedDEtaTRD%s",cutTM[i].Data()),
2922  Form("#Delta #eta of cluster-track vs cluster energy, SM behind TRD, %s",cutTM[i].Data()),
2923  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2924  fhTrackMatchedDEtaTRD[i]->SetYTitle("#Delta #eta");
2925  fhTrackMatchedDEtaTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2926 
2927  fhTrackMatchedDPhiTRD[i] = new TH2F
2928  (Form("hTrackMatchedDPhiTRD%s",cutTM[i].Data()),
2929  Form("#Delta #varphi of cluster-track vs cluster energy, SM behing TRD, %s",cutTM[i].Data()),
2930  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2931  fhTrackMatchedDPhiTRD[i]->SetYTitle("#Delta #varphi (rad)");
2932  fhTrackMatchedDPhiTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2933 
2934  fhEOverPTRD[i] = new TH2F
2935  (Form("hEOverPTRD%s",cutTM[i].Data()),
2936  Form("matched track #it{E}/#it{p} vs cluster E, behind TRD, %s",cutTM[i].Data()),
2937  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2938  fhEOverPTRD[i]->SetXTitle("#it{E} (GeV)");
2939  fhEOverPTRD[i]->SetYTitle("#it{E}/#it{p}");
2940 
2941  outputContainer->Add(fhTrackMatchedDEtaTRD[i]) ;
2942  outputContainer->Add(fhTrackMatchedDPhiTRD[i]) ;
2943  outputContainer->Add(fhEOverPTRD[i]);
2944  }
2945 
2946  if ( IsDataMC() )
2947  {
2949  (Form("hTrackMatchedDEtaMCNoOverlap%s",cutTM[i].Data()),
2950  Form("#Delta #eta of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
2951  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2952  fhTrackMatchedDEtaMCNoOverlap[i]->SetYTitle("#Delta #eta");
2953  fhTrackMatchedDEtaMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2954 
2956  (Form("hTrackMatchedDPhiMCNoOverlap%s",cutTM[i].Data()),
2957  Form("#Delta #varphi of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
2958  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2959  fhTrackMatchedDPhiMCNoOverlap[i]->SetYTitle("#Delta #varphi (rad)");
2960  fhTrackMatchedDPhiMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2961 
2962  outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[i]) ;
2963  outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[i]) ;
2965  (Form("hTrackMatchedDEtaMCOverlap%s",cutTM[i].Data()),
2966  Form("#Delta #eta of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
2967  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2968  fhTrackMatchedDEtaMCOverlap[i]->SetYTitle("#Delta #eta");
2969  fhTrackMatchedDEtaMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2970 
2972  (Form("hTrackMatchedDPhiMCOverlap%s",cutTM[i].Data()),
2973  Form("#Delta #varphi of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
2974  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2975  fhTrackMatchedDPhiMCOverlap[i]->SetYTitle("#Delta #varphi (rad)");
2976  fhTrackMatchedDPhiMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2977 
2978  outputContainer->Add(fhTrackMatchedDEtaMCOverlap[i]) ;
2979  outputContainer->Add(fhTrackMatchedDPhiMCOverlap[i]) ;
2980 
2982  (Form("hTrackMatchedDEtaMCConversion%s",cutTM[i].Data()),
2983  Form("#Delta #eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
2984  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2985  fhTrackMatchedDEtaMCConversion[i]->SetYTitle("#Delta #eta");
2986  fhTrackMatchedDEtaMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2987 
2989  (Form("hTrackMatchedDPhiMCConversion%s",cutTM[i].Data()),
2990  Form("#Delta #varphi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
2991  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2992  fhTrackMatchedDPhiMCConversion[i]->SetYTitle("#Delta #varphi (rad)");
2993  fhTrackMatchedDPhiMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2994 
2995  outputContainer->Add(fhTrackMatchedDEtaMCConversion[i]) ;
2996  outputContainer->Add(fhTrackMatchedDPhiMCConversion[i]) ;
2997 
2999  (Form("hTrackMatchedMCParticle%s_BeforeTM",cutTM[i].Data()),
3000  Form("Origin of particle vs energy %s, before TM criteria",cutTM[i].Data()),
3001  nptbins,ptmin,ptmax,11,0,11);
3002  fhTrackMatchedMCParticleBeforeTM[i]->SetXTitle("#it{E} (GeV)");
3003  //fhTrackMatchedMCParticleBeforeTM[i]->SetYTitle("Particle type");
3004  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3005  fhTrackMatchedMCParticleBeforeTM[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3006 
3007  outputContainer->Add(fhTrackMatchedMCParticleBeforeTM[i]);
3008 
3010  (Form("hTrackMatchedMCParticleTrackPt%s_BeforeTM",cutTM[i].Data()),
3011  Form("Origin of particle vs track #it{p}_{T} %s, before TM criteria",cutTM[i].Data()),
3012  nptbins,ptmin,ptmax,11,0,11);
3013  fhTrackMatchedMCParticleTrackPtBeforeTM[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
3014  //fhTrackMatchedMCParticleTrackPtBeforeTM[i]->SetYTitle("Particle type");
3015  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3016  fhTrackMatchedMCParticleTrackPtBeforeTM[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3017 
3018  outputContainer->Add(fhTrackMatchedMCParticleTrackPtBeforeTM[i]);
3019 
3021  (Form("hTrackMatchedMCParticle%s",cutTM[i].Data()),
3022  Form("Origin of particle vs energy %s",cutTM[i].Data()),
3023  nptbins,ptmin,ptmax,11,0,11);
3024  fhTrackMatchedMCParticle[i]->SetXTitle("#it{E} (GeV)");
3025  //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
3026  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3027  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3028 
3029  outputContainer->Add(fhTrackMatchedMCParticle[i]);
3030 
3032  (Form("hTrackMatchedMCParticleTrackPt%s",cutTM[i].Data()),
3033  Form("Origin of particle vs track #it{p}_{T} %s",cutTM[i].Data()),
3034  nptbins,ptmin,ptmax,11,0,11);
3035  fhTrackMatchedMCParticleTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
3036  //fhTrackMatchedMCParticleTrackPt[i]->SetYTitle("Particle type");
3037  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3038  fhTrackMatchedMCParticleTrackPt[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3039 
3040  outputContainer->Add(fhTrackMatchedMCParticleTrackPt[i]);
3041 
3043  (Form("hTrackMatchedMCParticle%s_Converted",cutTM[i].Data()),
3044  Form("Origin of particle vs energy %s, suffer conversion",cutTM[i].Data()),
3045  nptbins,ptmin,ptmax,11,0,11);
3046  fhTrackMatchedMCParticleConverted[i]->SetXTitle("#it{E} (GeV)");
3047  //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
3048  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3049  fhTrackMatchedMCParticleConverted[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3050 
3051  outputContainer->Add(fhTrackMatchedMCParticleConverted[i]);
3052 
3054  (Form("hTrackMatchedMCParticle%s_EOverP",cutTM[i].Data()),
3055  Form("Origin of particle vs energy %s vs #it{E}/#it{p}",cutTM[i].Data()),
3056  nptbins,ptmin,ptmax,11,0,11,100,0,20);
3057  fhTrackMatchedMCParticleVsEOverP[i]->SetXTitle("#it{E} (GeV)");
3058  fhTrackMatchedMCParticleVsEOverP[i]->SetZTitle("#it{E}/#it{p}");
3059  //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
3060  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3061  fhTrackMatchedMCParticleVsEOverP[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3062 
3063  outputContainer->Add(fhTrackMatchedMCParticleVsEOverP[i]);
3064 
3066  (Form("hTrackMatchedMCParticle%s_NOverlaps",cutTM[i].Data()),
3067  Form("Origin of particle vs energy %s vs number of particle overlaps",cutTM[i].Data()),
3068  nptbins,ptmin,ptmax,11,0,11,10,0,10);
3069  fhTrackMatchedMCParticleVsNOverlaps[i]->SetXTitle("#it{E} (GeV)");
3070  fhTrackMatchedMCParticleVsNOverlaps[i]->SetZTitle("#it{n}_{overlaps}");
3071  //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
3072  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3073  fhTrackMatchedMCParticleVsNOverlaps[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3074 
3075  outputContainer->Add(fhTrackMatchedMCParticleVsNOverlaps[i]);
3076 
3078  (Form("hTrackMatchedMCParticle%s_ErecEgenRat",cutTM[i].Data()),
3079  Form("Origin of particle vs energy %s vs #it{E}_{rec}/#it{E}_{gen}",cutTM[i].Data()),
3080  nptbins,ptmin,ptmax,11,0,11,100,0,2);
3081  fhTrackMatchedMCParticleVsErecEgen[i]->SetXTitle("#it{E}_{rec} (GeV)");
3082  fhTrackMatchedMCParticleVsErecEgen[i]->SetZTitle("#it{E}_{rec}/#it{E}_{gen}");
3083  //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
3084  for(Int_t imcpart = 1; imcpart < 12; imcpart++)
3085  fhTrackMatchedMCParticleVsErecEgen[i]->GetYaxis()->SetBinLabel(imcpart,mcPartLabels[imcpart-1]);
3086 
3087  outputContainer->Add(fhTrackMatchedMCParticleVsErecEgen[i]);
3088 
3089  }
3090  }
3091  }
3092 
3093  if(IsPileUpAnalysisOn())
3094  {
3095  TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
3096 
3097  for(Int_t i = 0 ; i < 7 ; i++)
3098  {
3099  fhPtPhotonPileUp[i] = new TH1F(Form("hPtPhotonPileUp%s",pileUpName[i].Data()),
3100  Form("Selected photon #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
3101  fhPtPhotonPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3102  outputContainer->Add(fhPtPhotonPileUp[i]);
3103 
3104  fhClusterTimeDiffPhotonPileUp[i] = new TH2F(Form("hClusterTimeDiffPhotonPileUp%s",pileUpName[i].Data()),
3105  Form("Photon cluster E vs #it{t}_{max}-#it{t}_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
3106  nptbins,ptmin,ptmax,400,-200,200);
3107  fhClusterTimeDiffPhotonPileUp[i]->SetXTitle("#it{E} (GeV)");
3108  fhClusterTimeDiffPhotonPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
3109  outputContainer->Add(fhClusterTimeDiffPhotonPileUp[i]);
3110  }
3111 
3112  fhTimePtPhotonNoCut = new TH2F ("hTimePtPhoton_NoCut","time of photon cluster vs pT of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3113  fhTimePtPhotonNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3114  fhTimePtPhotonNoCut->SetYTitle("#it{time} (ns)");
3115  outputContainer->Add(fhTimePtPhotonNoCut);
3116 
3117  fhTimePtPhotonSPD = new TH2F ("hTimePtPhoton_SPD","time of photon cluster vs pT of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3118  fhTimePtPhotonSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3119  fhTimePtPhotonSPD->SetYTitle("#it{time} (ns)");
3120  outputContainer->Add(fhTimePtPhotonSPD);
3121 
3122  fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,20,0,20);
3123  fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
3124  fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
3125  outputContainer->Add(fhTimeNPileUpVertSPD);
3126 
3127  fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 20,0,20 );
3128  fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
3129  fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
3130  outputContainer->Add(fhTimeNPileUpVertTrack);
3131 
3132  fhPtPhotonNPileUpSPDVtx = new TH2F ("hPtPhoton_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
3133  nptbins,ptmin,ptmax,20,0,20);
3134  fhPtPhotonNPileUpSPDVtx->SetYTitle("# vertex ");
3135  fhPtPhotonNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3136  outputContainer->Add(fhPtPhotonNPileUpSPDVtx);
3137 
3138  fhPtPhotonNPileUpTrkVtx = new TH2F ("hPtPhoton_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
3139  nptbins,ptmin,ptmax, 20,0,20 );
3140  fhPtPhotonNPileUpTrkVtx->SetYTitle("# vertex ");
3141  fhPtPhotonNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3142  outputContainer->Add(fhPtPhotonNPileUpTrkVtx);
3143 
3144  fhPtPhotonNPileUpSPDVtxTimeCut = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
3145  nptbins,ptmin,ptmax,20,0,20);
3146  fhPtPhotonNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
3147  fhPtPhotonNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3148  outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut);
3149 
3150  fhPtPhotonNPileUpTrkVtxTimeCut = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
3151  nptbins,ptmin,ptmax, 20,0,20 );
3152  fhPtPhotonNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
3153  fhPtPhotonNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3154  outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut);
3155 
3156  fhPtPhotonNPileUpSPDVtxTimeCut2 = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
3157  nptbins,ptmin,ptmax,20,0,20);
3158  fhPtPhotonNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
3159  fhPtPhotonNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3160  outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut2);
3161 
3162  fhPtPhotonNPileUpTrkVtxTimeCut2 = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
3163  nptbins,ptmin,ptmax, 20,0,20 );
3164  fhPtPhotonNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
3165  fhPtPhotonNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3166  outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut2);
3167 
3168  }
3169 
3171  {
3172  for(Int_t ieta = 0; ieta < 4; ieta++)
3173  {
3174  for(Int_t iphi = 0; iphi < 3; iphi++)
3175  {
3176 // fhLam0EMCALRegion[ieta][iphi] =
3177 // new TH2F(Form("hLam0_eta%d_phi%d",ieta,iphi),
3178 // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, region eta %d, phi %d",ieta,iphi),
3179 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3180 // fhLam0EMCALRegion[ieta][iphi]->SetYTitle("#lambda_{0}^{2}");
3181 // fhLam0EMCALRegion[ieta][iphi]->SetXTitle("#it{p}_{T} (GeV)");
3182 // outputContainer->Add(fhLam0EMCALRegion[ieta][iphi]) ;
3183 //
3184 // if(GetFirstSMCoveredByTRD() >= 0)
3185 // {
3186 // fhLam0EMCALRegionTRD[ieta][iphi] =
3187 // new TH2F(Form("hLam0TRD_eta%d_phi%d",ieta,iphi),
3188 // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, region eta %d, phi %d, SM covered by TRD",ieta,iphi),
3189 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3190 // fhLam0EMCALRegionTRD[ieta][iphi]->SetYTitle("#lambda_{0}^{2}");
3191 // fhLam0EMCALRegionTRD[ieta][iphi]->SetXTitle("#it{p}_{T} (GeV)");
3192 // outputContainer->Add(fhLam0EMCALRegionTRD[ieta][iphi]) ;
3193 // } // TRD
3194 
3195  for(Int_t ism = 0; ism < fNModules; ism++)
3196  {
3197  if(ism < fFirstModule || ism > fLastModule) continue;
3198 
3199  fhLam0EMCALRegionPerSM[ieta][iphi][ism] =
3200  new TH2F(Form("hLam0_eta%d_phi%d_sm%d",ieta,iphi,ism),
3201  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, sm %d, region eta %d, phi %d",ism,ieta,iphi),
3202  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3203  fhLam0EMCALRegionPerSM[ieta][iphi][ism]->SetYTitle("#lambda_{0}^{2}");
3204  fhLam0EMCALRegionPerSM[ieta][iphi][ism]->SetXTitle("#it{p}_{T} (GeV)");
3205  outputContainer->Add(fhLam0EMCALRegionPerSM[ieta][iphi][ism]) ;
3206 
3207  fhLam1EMCALRegionPerSM[ieta][iphi][ism] =
3208  new TH2F(Form("hLam1_eta%d_phi%d_sm%d",ieta,iphi,ism),
3209  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{1}^{2}, sm %d, region eta %d, phi %d",ism,ieta,iphi),
3210  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3211  fhLam1EMCALRegionPerSM[ieta][iphi][ism]->SetYTitle("#lambda_{1}^{2}");
3212  fhLam1EMCALRegionPerSM[ieta][iphi][ism]->SetXTitle("#it{p}_{T} (GeV)");
3213  outputContainer->Add(fhLam1EMCALRegionPerSM[ieta][iphi][ism]) ;
3214  }
3215  } // iphi
3216  } // ieta
3217 
3218  Float_t ptLimit[] = {2,3,4,5,6,8,10,12};
3219  TString l0bin [] = {"0.23<#lambda^{2}_{0}<0.26","0.3<#lambda^{2}_{0}<0.4"};
3220 
3221  for(Int_t il0 = 0; il0 < 2; il0++)
3222  {
3223  for(Int_t ipt = 0; ipt < 7; ipt++)
3224  {
3225  fhEtaPhiLam0BinPtBin[il0][ipt] = new TH2F
3226  (Form("hEtaPhiLam0Bin%d_PtBin%d",il0,ipt),
3227  Form("#eta vs #varphi in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s",
3228  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3229  netabins,etamin,etamax,nphibins,phimin,phimax);
3230  fhEtaPhiLam0BinPtBin[il0][ipt]->SetYTitle("#varphi (rad)");
3231  fhEtaPhiLam0BinPtBin[il0][ipt]->SetXTitle("#eta");
3232  outputContainer->Add(fhEtaPhiLam0BinPtBin[il0][ipt]) ;
3233 
3234  fhColRowLam0BinPtBin[il0][ipt] = new TH2F
3235  (Form("hColRowLam0Bin%d_PtBin%d",il0,ipt),
3236  Form("row vs column in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s, w > 0",
3237  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3238  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3239  fhColRowLam0BinPtBin[il0][ipt]->SetYTitle("row");
3240  fhColRowLam0BinPtBin[il0][ipt]->SetXTitle("column");
3241  outputContainer->Add(fhColRowLam0BinPtBin[il0][ipt]) ;
3242 
3243  fhColRowLam0BinPtBinWeighted[il0][ipt] = new TH2F
3244  (Form("hColRowLam0Bin%d_PtBin%dWeighted",il0,ipt),
3245  Form("cluster cell row vs column weighted in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s",
3246  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3247  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3248  fhColRowLam0BinPtBinWeighted[il0][ipt]->SetYTitle("row");
3249  fhColRowLam0BinPtBinWeighted[il0][ipt]->SetXTitle("column");
3250  outputContainer->Add(fhColRowLam0BinPtBinWeighted[il0][ipt]) ;
3251 
3252  fhColRowLam0BinPtBinLargeTime[il0][ipt] = new TH2F
3253  (Form("hColRowLam0Bin%d_PtBin%d_LargeTimeInClusterCell",il0,ipt),
3254  Form("row vs column in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s, |t| > 50 ns, w > 0",
3255  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3256  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3257  fhColRowLam0BinPtBinLargeTime[il0][ipt]->SetYTitle("row");
3258  fhColRowLam0BinPtBinLargeTime[il0][ipt]->SetXTitle("column");
3259  outputContainer->Add(fhColRowLam0BinPtBinLargeTime[il0][ipt]) ;
3260 
3261 // fhEtaPhiLam0BinPtBinSMShared[il0][ipt] = new TH2F
3262 // (Form("hEtaPhiLam0Bin%d_PtBin%d_SMShared",il0,ipt),
3263 // Form("#eta vs #varphi in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s",
3264 // ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3265 // netabins,etamin,etamax,nphibins,phimin,phimax);
3266 // fhEtaPhiLam0BinPtBinSMShared[il0][ipt]->SetYTitle("#varphi (rad)");
3267 // fhEtaPhiLam0BinPtBinSMShared[il0][ipt]->SetXTitle("#eta");
3268 // outputContainer->Add(fhEtaPhiLam0BinPtBinSMShared[il0][ipt]) ;
3269 //
3270 // fhColRowLam0BinPtBinSMShared[il0][ipt] = new TH2F
3271 // (Form("hColRowLam0Bin%d_PtBin%d_SMShared",il0,ipt),
3272 // Form("row vs column in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s, w > 0",
3273 // ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3274 // ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3275 // fhColRowLam0BinPtBinSMShared[il0][ipt]->SetYTitle("row");
3276 // fhColRowLam0BinPtBinSMShared[il0][ipt]->SetXTitle("column");
3277 // outputContainer->Add(fhColRowLam0BinPtBinSMShared[il0][ipt]) ;
3278 
3279  fhEtaPhiLargeTimeInClusterCell[il0][ipt] = new TH2F
3280  (Form("hEtaPhiLam0Bin%d_PtBin%d_LargeTimeInClusterCell",il0,ipt),
3281  Form("#eta vs #varphi in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, t > 50 ns, %s",
3282  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3283  netabins,etamin,etamax,nphibins,phimin,phimax);
3284  fhEtaPhiLargeTimeInClusterCell[il0][ipt]->SetYTitle("#varphi (rad)");
3285  fhEtaPhiLargeTimeInClusterCell[il0][ipt]->SetXTitle("#eta");
3286  outputContainer->Add(fhEtaPhiLargeTimeInClusterCell[il0][ipt]) ;
3287 
3288  fhCellClusterIndexEAndTime[il0][ipt] = new TH2F
3289  (Form("hCellClusterIndexEAndTimeLam0Bin%d_PtBin%d",il0,ipt),
3290  Form("#it{t}_{cell} vs cell index (E sorted) in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, %s",
3291  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3292  ntimebins,timemin,timemax,30,0,30);
3293  fhCellClusterIndexEAndTime[il0][ipt] ->SetYTitle("Index (sorted #it{E})");
3294  fhCellClusterIndexEAndTime[il0][ipt] ->SetXTitle("#it{time}_{cell} (ns)");
3295  outputContainer->Add(fhCellClusterIndexEAndTime[il0][ipt] ) ;
3296 
3297  fhCellClusterEAndTime[il0][ipt] = new TH2F
3298  (Form("hCellClusterEAndTimeLam0Bin%d_PtBin%d",il0,ipt),
3299  Form("#it{E}_{cell} vs #it{t}_{cell} in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, %s",
3300  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3301  ntimebins,timemin,timemax,100,0,5);
3302  fhCellClusterEAndTime[il0][ipt] ->SetYTitle("#it{E}_{cell} (GeV)");
3303  fhCellClusterEAndTime[il0][ipt] ->SetXTitle("#it{time}_{cell} (ns)");
3304  outputContainer->Add(fhCellClusterEAndTime[il0][ipt] ) ;
3305 
3306  fhCellClusterEFracAndTime[il0][ipt] = new TH2F
3307  (Form("hCellClusterEFracAndTimeLam0Bin%d_PtBin%d",il0,ipt),
3308  Form("#it{E}_{cell} vs #it{t}_{cell} in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, %s",
3309  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3310  ntimebins,timemin,timemax,100,0,1);
3311  fhCellClusterEFracAndTime[il0][ipt] ->SetYTitle("#it{E}_{cell} (GeV)");
3312  fhCellClusterEFracAndTime[il0][ipt] ->SetXTitle("#it{time}_{cell} (ns)");
3313  outputContainer->Add(fhCellClusterEFracAndTime[il0][ipt] ) ;
3314  } // pt bin
3315 
3316  for(Int_t ism = 0; ism < fNModules; ism++)
3317  {
3318  if(ism < fFirstModule || ism > fLastModule) continue;
3319 
3320  fhLam1Lam0BinPerSM[il0][ism] = new TH2F
3321  (Form("hLam1Lam0Bin%d_sm%d",il0,ism),
3322  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, %s",ism,l0bin[il0].Data()),
3323  nptbins,ptmin,ptmax,40,0,0.4);
3324  fhLam1Lam0BinPerSM[il0][ism]->SetYTitle("#lambda^{2}_{1}");
3325  fhLam1Lam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3326  outputContainer->Add(fhLam1Lam0BinPerSM[il0][ism]) ;
3327 
3328  fhTimeLam0BinPerSM[il0][ism] = new TH2F
3329  (Form("hTimeLam0Bin%d_sm%d",il0,ism),
3330  Form("#it{p}_{T} vs cluster cell time in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3331  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3332  fhTimeLam0BinPerSM[il0][ism]->SetYTitle("cell time (ns)");
3333  fhTimeLam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3334  outputContainer->Add(fhTimeLam0BinPerSM[il0][ism]) ;
3335 
3336  fhTimeLam0BinPerSMWeighted[il0][ism] = new TH2F
3337  (Form("hTimeLam0Bin%d_sm%d_Weighted",il0,ism),
3338  Form("#it{p}_{T} vs cluster cell time weighted in sm %d, %s",ism,l0bin[il0].Data()),
3339  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3340  fhTimeLam0BinPerSMWeighted[il0][ism]->SetYTitle("cell time (ns)");
3341  fhTimeLam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3342  outputContainer->Add(fhTimeLam0BinPerSMWeighted[il0][ism]) ;
3343 
3344  fhDTimeLam0BinPerSM[il0][ism] = new TH2F
3345  (Form("hDTimeLam0Bin%d_sm%d",il0,ism),
3346  Form("#it{p}_{T} vs t_{cluster}-t_{cell} in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3347  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3348  fhDTimeLam0BinPerSM[il0][ism]->SetYTitle("t_{cluster}-t_{cell} (ns)");
3349  fhDTimeLam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3350  outputContainer->Add(fhDTimeLam0BinPerSM[il0][ism]) ;
3351 
3352  fhDTimeLam0BinPerSMWeighted[il0][ism] = new TH2F
3353  (Form("hDTimeLam0Bin%d_sm%d_Weighted",il0,ism),
3354  Form("#it{p}_{T} vs t_{cluster}-t_{cell} weighted in sm %d, %s",ism,l0bin[il0].Data()),
3355  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3356  fhDTimeLam0BinPerSMWeighted[il0][ism]->SetYTitle("t_{cluster}-t_{cell} (ns)");
3357  fhDTimeLam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3358  outputContainer->Add(fhDTimeLam0BinPerSMWeighted[il0][ism]) ;
3359 
3360  fhCellClusterEFracLam0BinPerSM[il0][ism] = new TH2F
3361  (Form("hCellClusterEFracLam0Bin%d_sm%d",il0,ism),
3362  Form("#it{p}_{T} vs cell E / cluster E in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3363  nptbins,ptmin,ptmax,100,0,1);
3364  fhCellClusterEFracLam0BinPerSM[il0][ism]->SetYTitle("cell E / cluster E ");
3365  fhCellClusterEFracLam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3366  outputContainer->Add(fhCellClusterEFracLam0BinPerSM[il0][ism]) ;
3367 
3368 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism] = new TH2F
3369 // (Form("hCellClusterEFracLam0Bin%d_sm%d_Weighted",il0,ism),
3370 // Form("#it{p}_{T} vs cell E / cluster E weighted in sm %d, %s",ism,l0bin[il0].Data()),
3371 // nptbins,ptmin,ptmax,100,0,1);
3372 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism]->SetYTitle("cell E / cluster E ");
3373 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3374 // outputContainer->Add(fhCellClusterEFracLam0BinPerSMWeighted[il0][ism]) ;
3375 
3377  (Form("hCellClusterEFracLam0Bin%d_sm%d_LargeTimeInClusterCell",il0,ism),
3378  Form("#it{p}_{T} vs cell E / cluster E in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3379  nptbins,ptmin,ptmax,100,0,1);
3380  fhCellClusterEFracLam0BinPerSMLargeTime[il0][ism]->SetYTitle("cell E / cluster E ");
3381  fhCellClusterEFracLam0BinPerSMLargeTime[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3382  outputContainer->Add(fhCellClusterEFracLam0BinPerSMLargeTime[il0][ism]) ;
3383 
3385  (Form("hCellClusterEFracLam0Bin%d_sm%d_LargeTimeInClusterCell_Total",il0,ism),
3386  Form("#it{p}_{T} vs cell E / cluster E in sm %d, %s, w > 0, all |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3387  nptbins,ptmin,ptmax,100,0,1);
3388  fhCellClusterEFracLam0BinPerSMLargeTimeTotal[il0][ism]->SetYTitle("cell E / cluster E ");
3389  fhCellClusterEFracLam0BinPerSMLargeTimeTotal[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3390  outputContainer->Add(fhCellClusterEFracLam0BinPerSMLargeTimeTotal[il0][ism]) ;
3391 
3392  fhCellClusterELam0BinPerSM[il0][ism] = new TH2F
3393  (Form("hCellClusterELam0Bin%d_sm%d",il0,ism),
3394  Form("#it{p}_{T} vs cell E in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3395  nptbins,ptmin,ptmax,500,0,10);
3396  fhCellClusterELam0BinPerSM[il0][ism]->SetYTitle("cell E (GeV)");
3397  fhCellClusterELam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3398  outputContainer->Add(fhCellClusterELam0BinPerSM[il0][ism]) ;
3399 
3401  (Form("hCellClusterELam0Bin%d_sm%d_Weighted",il0,ism),
3402  Form("#it{p}_{T} vs cell E weighted in sm %d, %s",ism,l0bin[il0].Data()),
3403  nptbins,ptmin,ptmax,500,0,10);
3404  fhCellClusterELam0BinPerSMWeighted[il0][ism]->SetYTitle("cell E (GeV)");
3405  fhCellClusterELam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3406  outputContainer->Add(fhCellClusterELam0BinPerSMWeighted[il0][ism]) ;
3407 
3409  (Form("hCellClusterELam0Bin%d_sm%d_LargeTimeInClusterCell",il0,ism),
3410  Form("#it{p}_{T} vs cell E in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3411  nptbins,ptmin,ptmax,500,0,10);
3412  fhCellClusterELam0BinPerSMLargeTime[il0][ism]->SetYTitle("cell E (GeV)");
3413  fhCellClusterELam0BinPerSMLargeTime[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3414  outputContainer->Add(fhCellClusterELam0BinPerSMLargeTime[il0][ism]) ;
3415 
3417  (Form("hCellClusterIndexELam0Bin%d_sm%d_LargeTimeInClusterCell",il0,ism),
3418  Form("#it{p}_{T} vs cell index (E sorted) in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3419  nptbins,ptmin,ptmax,30,0,30);
3420  fhCellClusterIndexELam0BinPerSMLargeTime[il0][ism]->SetYTitle("Index");
3421  fhCellClusterIndexELam0BinPerSMLargeTime[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3422  outputContainer->Add(fhCellClusterIndexELam0BinPerSMLargeTime[il0][ism]) ;
3423 
3424  fhNCellsWithLargeTimeInCluster[il0][ism] = new TH2F
3425  (Form("hNCellsWithLargeTimeInClusterLam0Bin%d_sm%d",il0,ism),
3426  Form("#it{p}_{T} vs number of cells in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3427  nptbins,ptmin,ptmax,30,0,30);
3428  fhNCellsWithLargeTimeInCluster[il0][ism]->SetYTitle("#it{n}_{cells}");
3429  fhNCellsWithLargeTimeInCluster[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3430  outputContainer->Add(fhNCellsWithLargeTimeInCluster[il0][ism]) ;
3431 
3432 // if(ism < 12)
3433 // {
3434 // fhTimeLam0BinPerSMShared[il0][ism] = new TH2F
3435 // (Form("hTimeLam0Bin%d_sm%d_SMShared",il0,ism),
3436 // Form("#it{p}_{T} vs cluster cell time in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3437 // nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3438 // fhTimeLam0BinPerSMShared[il0][ism]->SetYTitle("cell time (ns)");
3439 // fhTimeLam0BinPerSMShared[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3440 // outputContainer->Add(fhTimeLam0BinPerSMShared[il0][ism]) ;
3441 // } // Run 1 SM
3442  } // sm
3443  } // l0 bin
3444 
3445  for(Int_t ism = 0; ism < fNModules; ism++)
3446  {
3447  if(ism < fFirstModule || ism > fLastModule) continue;
3448 
3450  (Form("hLam0_sm%d_LargeTimeInClusterCell",ism),
3451  Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d,|t_{secondary cell}| > 50 ns",ism),
3452  nptbins,ptmin,ptmax,40,0,0.4);
3453  fhLam0PerSMLargeTimeInClusterCell[ism]->SetYTitle("#lambda^{2}_{0}");
3454  fhLam0PerSMLargeTimeInClusterCell[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3455  outputContainer->Add(fhLam0PerSMLargeTimeInClusterCell[ism]) ;
3456 
3458  (Form("hLam1_sm%d_LargeTimeInClusterCell",ism),
3459  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, |t_{secondary cell}| > 50 ns",ism),
3460  nptbins,ptmin,ptmax,40,0,0.4);
3461  fhLam1PerSMLargeTimeInClusterCell[ism]->SetYTitle("#lambda^{2}_{1}");
3462  fhLam1PerSMLargeTimeInClusterCell[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3463  outputContainer->Add(fhLam1PerSMLargeTimeInClusterCell[ism]) ;
3464 
3465 // fhLam0PerSMSPDPileUp[ism] = new TH2F
3466 // (Form("hLam0_sm%d_SPDPileUp",ism),
3467 // Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d, Pile-up event SPD",ism),
3468 // nptbins,ptmin,ptmax,40,0,0.4);
3469 // fhLam0PerSMSPDPileUp[ism]->SetYTitle("#lambda^{2}_{0}");
3470 // fhLam0PerSMSPDPileUp[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3471 // outputContainer->Add(fhLam0PerSMSPDPileUp[ism]) ;
3472 //
3473 // fhLam1PerSMSPDPileUp[ism] = new TH2F
3474 // (Form("hLam1_sm%d_SPDPileUp",ism),
3475 // Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, Pile-up event SPD",ism),
3476 // nptbins,ptmin,ptmax,40,0,0.4);
3477 // fhLam1PerSMSPDPileUp[ism]->SetYTitle("#lambda^{2}_{1}");
3478 // fhLam1PerSMSPDPileUp[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3479 // outputContainer->Add(fhLam1PerSMSPDPileUp[ism]) ;
3480 
3481 // if(ism < 12)
3482 // {
3483 // fhLam0PerSMShared[ism] = new TH2F
3484 // (Form("hLam0_sm%d_SMShared",ism),
3485 // Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d, SM shared",ism),
3486 // nptbins,ptmin,ptmax,40,0,0.4);
3487 // fhLam0PerSMShared[ism]->SetYTitle("#lambda^{2}_{0}");
3488 // fhLam0PerSMShared[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3489 // outputContainer->Add(fhLam0PerSMShared[ism]) ;
3490 //
3491 // fhLam1PerSMShared[ism] = new TH2F
3492 // (Form("hLam1_sm%d_SMShared",ism),
3493 // Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, SM shared",ism),
3494 // nptbins,ptmin,ptmax,40,0,0.4);
3495 // fhLam1PerSMShared[ism]->SetYTitle("#lambda^{2}_{1}");
3496 // fhLam1PerSMShared[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3497 // outputContainer->Add(fhLam1PerSMShared[ism]) ;
3498 // } // run1
3499  } // sm
3500 
3501  for(Int_t ilarge = 0; ilarge < 5; ilarge++)
3502  {
3504  (Form("hLam0_NLargeTimeInClusterCell%d",ilarge),
3505  Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d,|t_{secondary cell}| > 50 ns",ilarge),
3506  nptbins,ptmin,ptmax,40,0,0.4);
3507  fhLam0PerNLargeTimeInClusterCell[ilarge]->SetYTitle("#lambda^{2}_{0}");
3508  fhLam0PerNLargeTimeInClusterCell[ilarge]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3509  outputContainer->Add(fhLam0PerNLargeTimeInClusterCell[ilarge]) ;
3510 
3512  (Form("hLam1_NLargeTimeInClusterCell%d",ilarge),
3513  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, |t_{secondary cell}| > 50 ns",ilarge),
3514  nptbins,ptmin,ptmax,40,0,0.4);
3515  fhLam1PerNLargeTimeInClusterCell[ilarge]->SetYTitle("#lambda^{2}_{1}");
3516  fhLam1PerNLargeTimeInClusterCell[ilarge]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3517  outputContainer->Add(fhLam1PerNLargeTimeInClusterCell[ilarge]) ;
3518  }
3519 
3520  } // regions in EMCal
3521 
3522 
3523  if ( IsDataMC() )
3524  {
3525  TString ptype[] = { "#gamma" , "#gamma_{#pi decay}" , "#gamma_{#eta decay}", "#gamma_{other decay}",
3526  "#pi^{0}" , "#eta" , "e^{#pm}" , "#gamma->e^{#pm}" ,
3527  "hadron?" , "Anti-N" , "Anti-P" ,
3528  "Neutron" , "Proton" , "#pi^{#pm}" ,
3529  "#gamma_{prompt}", "#gamma_{fragmentation}", "#gamma_{ISR}" , "String" } ;
3530 
3531  TString pname[] = { "Photon" , "PhotonPi0Decay" , "PhotonEtaDecay", "PhotonOtherDecay",
3532  "Pi0" , "Eta" , "Electron" , "Conversion" ,
3533  "Hadron" , "AntiNeutron" , "AntiProton" ,
3534  "Neutron" , "Proton" , "ChPion" ,
3535  "PhotonPrompt", "PhotonFragmentation", "PhotonISR" , "String" } ;
3536 
3537  for(Int_t i = 0; i < fNOriginHistograms; i++)
3538  {
3539  fhMCE[i] = new TH1F(Form("hE_MC%s",pname[i].Data()),
3540  Form("cluster from %s : E ",ptype[i].Data()),
3541  nptbins,ptmin,ptmax);
3542  fhMCE[i]->SetXTitle("#it{E} (GeV)");
3543  outputContainer->Add(fhMCE[i]) ;
3544 
3545  fhMCPt[i] = new TH1F(Form("hPt_MC%s",pname[i].Data()),
3546  Form("cluster from %s : #it{p}_{T} ",ptype[i].Data()),
3547  nptbins,ptmin,ptmax);
3548  fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3549  outputContainer->Add(fhMCPt[i]) ;
3550 
3551  fhMCEta[i] = new TH2F(Form("hEta_MC%s",pname[i].Data()),
3552  Form("cluster from %s : #eta ",ptype[i].Data()),
3553  nptbins,ptmin,ptmax,netabins,etamin,etamax);
3554  fhMCEta[i]->SetYTitle("#eta");
3555  fhMCEta[i]->SetXTitle("#it{E} (GeV)");
3556  outputContainer->Add(fhMCEta[i]) ;
3557 
3558  fhMCPhi[i] = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
3559  Form("cluster from %s : #varphi ",ptype[i].Data()),
3560  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3561  fhMCPhi[i]->SetYTitle("#varphi (rad)");
3562  fhMCPhi[i]->SetXTitle("#it{E} (GeV)");
3563  outputContainer->Add(fhMCPhi[i]) ;
3564 
3565 
3566  fhMCDeltaE[i] = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
3567  Form("MC - Reco E from %s",pname[i].Data()),
3568  nptbins,ptmin,ptmax, 200,-50,50);
3569  fhMCDeltaE[i]->SetYTitle("#Delta #it{E} (GeV)");
3570  fhMCDeltaE[i]->SetXTitle("#it{E} (GeV)");
3571  outputContainer->Add(fhMCDeltaE[i]);
3572 
3573  fhMCDeltaPt[i] = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
3574  Form("MC - Reco #it{p}_{T} from %s",pname[i].Data()),
3575  nptbins,ptmin,ptmax, 200,-50,50);
3576  fhMCDeltaPt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
3577  fhMCDeltaPt[i]->SetYTitle("#Delta #it{p}_{T} (GeV/#it{c})");
3578  outputContainer->Add(fhMCDeltaPt[i]);
3579 
3580  fhMC2E[i] = new TH2F (Form("h2E_MC%s",pname[i].Data()),
3581  Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
3582  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3583  fhMC2E[i]->SetXTitle("#it{E}_{rec} (GeV)");
3584  fhMC2E[i]->SetYTitle("#it{E}_{gen} (GeV)");
3585  outputContainer->Add(fhMC2E[i]);
3586 
3587  fhMC2Pt[i] = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
3588  Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
3589  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3590  fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
3591  fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/#it{c})");
3592  outputContainer->Add(fhMC2Pt[i]);
3593  }
3594 
3595  TString pptype[] = { "#gamma" , "#gamma_{#pi decay}" ,
3596  "#gamma_{#eta decay}", "#gamma_{other decay}" ,
3597  "#gamma_{prompt}" , "#gamma_{fragmentation}", "#gamma_{ISR}" } ;
3598 
3599  TString ppname[] = { "Photon" , "PhotonPi0Decay" ,
3600  "PhotonEtaDecay", "PhotonOtherDecay" ,
3601  "PhotonPrompt" , "PhotonFragmentation", "PhotonISR" } ;
3602 
3603  for(Int_t i = 0; i < fNPrimaryHistograms; i++)
3604  {
3605  fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
3606  Form("primary photon %s : E ",pptype[i].Data()),
3607  nptbins,ptmin,ptmax);
3608  fhEPrimMC[i]->SetXTitle("#it{E} (GeV)");
3609  outputContainer->Add(fhEPrimMC[i]) ;
3610 
3611  fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
3612  Form("primary photon %s : #it{p}_{T} ",pptype[i].Data()),
3613  nptbins,ptmin,ptmax);
3614  fhPtPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3615  outputContainer->Add(fhPtPrimMC[i]) ;
3616 
3617  fhYPrimMC[i] = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
3618  Form("primary photon %s : Rapidity ",pptype[i].Data()),
3619  nptbins,ptmin,ptmax,200,-2,2);
3620  fhYPrimMC[i]->SetYTitle("Rapidity");
3621  fhYPrimMC[i]->SetXTitle("#it{E} (GeV)");
3622  outputContainer->Add(fhYPrimMC[i]) ;
3623 
3624  fhEtaPrimMC[i] = new TH2F(Form("hEtaPrim_MC%s",ppname[i].Data()),
3625  Form("primary photon %s : #eta",pptype[i].Data()),
3626  nptbins,ptmin,ptmax,200,-2,2);
3627  fhEtaPrimMC[i]->SetYTitle("#eta");
3628  fhEtaPrimMC[i]->SetXTitle("#it{E} (GeV)");
3629  outputContainer->Add(fhEtaPrimMC[i]) ;
3630 
3631  fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
3632  Form("primary photon %s : #varphi ",pptype[i].Data()),
3633  nptbins,ptmin,ptmax,nphibins,0,TMath::TwoPi());
3634  fhPhiPrimMC[i]->SetYTitle("#varphi (rad)");
3635  fhPhiPrimMC[i]->SetXTitle("#it{E} (GeV)");
3636  outputContainer->Add(fhPhiPrimMC[i]) ;
3637 
3638 
3639  fhEPrimMCAcc[i] = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
3640  Form("primary photon %s in acceptance: E ",pptype[i].Data()),
3641  nptbins,ptmin,ptmax);
3642  fhEPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3643  outputContainer->Add(fhEPrimMCAcc[i]) ;
3644 
3645  fhPtPrimMCAcc[i] = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
3646  Form("primary photon %s in acceptance: #it{p}_{T} ",pptype[i].Data()),
3647  nptbins,ptmin,ptmax);
3648  fhPtPrimMCAcc[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3649  outputContainer->Add(fhPtPrimMCAcc[i]) ;
3650 
3651  fhYPrimMCAcc[i] = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
3652  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
3653  nptbins,ptmin,ptmax,100,-1,1);
3654  fhYPrimMCAcc[i]->SetYTitle("Rapidity");
3655  fhYPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3656  outputContainer->Add(fhYPrimMCAcc[i]) ;
3657 
3658  fhEtaPrimMCAcc[i] = new TH2F(Form("hEtaPrimAcc_MC%s",ppname[i].Data()),
3659  Form("primary photon %s in acceptance: #eta ",pptype[i].Data()),
3660  nptbins,ptmin,ptmax,netabins,etamin,etamax);
3661  fhEtaPrimMCAcc[i]->SetYTitle("#eta");
3662  fhEtaPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3663  outputContainer->Add(fhEtaPrimMCAcc[i]) ;
3664 
3665  fhPhiPrimMCAcc[i] = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
3666  Form("primary photon %s in acceptance: #varphi ",pptype[i].Data()),
3667  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3668  fhPhiPrimMCAcc[i]->SetYTitle("#varphi (rad)");
3669  fhPhiPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3670  outputContainer->Add(fhPhiPrimMCAcc[i]) ;
3671  }
3672 
3673  if ( fFillSSHistograms )
3674  {
3675  TString ptypess[] = { "#gamma","#gamma->e^{#pm}",
3676  "#pi^{0}","#pi^{0}->#gamma->e^{#pm}",
3677  "#eta" ,"#eta->#gamma->e^{#pm}"
3678  ,"e^{#pm}","other"} ;
3679 
3680  TString pnamess[] = { "Photon","PhotonConv",
3681  "Pi0","Pi0Conv","Eta","EtaConv","Electron","Other"} ;
3682 
3683  for(Int_t i = 0; i < fgkNssTypes; i++)
3684  {
3685  fhMCELambda0[i] = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
3686  Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
3687  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3688  fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
3689  fhMCELambda0[i]->SetXTitle("#it{E} (GeV)");
3690  outputContainer->Add(fhMCELambda0[i]) ;
3691 
3692  fhMCPtLambda0[i] = new TH2F(Form("hPtLambda0_MC%s",pnamess[i].Data()),
3693  Form("cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptypess[i].Data()),
3694  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3695  fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
3696  fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3697  outputContainer->Add(fhMCPtLambda0[i]) ;
3698 
3699  if(!GetReader()->IsEmbeddedClusterSelectionOn())
3700  {
3701  for(Int_t iover = 0; iover < 3; iover++)
3702  {
3703  fhMCPtLambda0Overlaps[i][iover] = new TH2F(Form("hPtLambda0_MC%s_Overlap%d",pnamess[i].Data(),iover),
3704  Form("cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, N Overlaps = %d",ptypess[i].Data(),iover),
3705  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3706  fhMCPtLambda0Overlaps[i][iover]->SetYTitle("#lambda_{0}^{2}");
3707  fhMCPtLambda0Overlaps[i][iover]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3708  outputContainer->Add(fhMCPtLambda0Overlaps[i][iover]) ;
3709  }
3710  }
3711 
3712  fhMCELambda1[i] = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
3713  Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
3714  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3715  fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
3716  fhMCELambda1[i]->SetXTitle("#it{E} (GeV)");
3717  outputContainer->Add(fhMCELambda1[i]) ;
3718 
3719  fhMCNCellsE[i] = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
3720  Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
3721  nptbins,ptmin,ptmax, nbins,nmin,nmax);
3722  fhMCNCellsE[i]->SetXTitle("#it{E} (GeV)");
3723  fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
3724  outputContainer->Add(fhMCNCellsE[i]);
3725 
3727  {
3728  fhMCEDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
3729  Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
3730  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3731  fhMCEDispersion[i]->SetYTitle("D^{2}");
3732  fhMCEDispersion[i]->SetXTitle("#it{E} (GeV)");
3733  outputContainer->Add(fhMCEDispersion[i]) ;
3734 
3735  fhMCMaxCellDiffClusterE[i] = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
3736  Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
3737  nptbins,ptmin,ptmax, 500,0,1.);
3738  fhMCMaxCellDiffClusterE[i]->SetXTitle("#it{E}_{cluster} (GeV) ");
3739  fhMCMaxCellDiffClusterE[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3740  outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
3741 
3742  fhMCLambda0vsClusterMaxCellDiffE0[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
3743  Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
3744  ssbins,ssmin,ssmax,500,0,1.);
3745  fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
3746  fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3747  outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
3748 
3749  fhMCLambda0vsClusterMaxCellDiffE2[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
3750  Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
3751  ssbins,ssmin,ssmax,500,0,1.);
3752  fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
3753  fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3754  outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
3755 
3756  fhMCLambda0vsClusterMaxCellDiffE6[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
3757  Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
3758  ssbins,ssmin,ssmax,500,0,1.);
3759  fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
3760  fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3761  outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
3762 
3763  fhMCNCellsvsClusterMaxCellDiffE0[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
3764  Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
3765  nbins/5,nmin,nmax/5,500,0,1.);
3766  fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
3767  fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3768  outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
3769 
3770  fhMCNCellsvsClusterMaxCellDiffE2[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
3771  Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
3772  nbins/5,nmin,nmax/5,500,0,1.);
3773  fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
3774  fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3775  outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
3776 
3777  fhMCNCellsvsClusterMaxCellDiffE6[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
3778  Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
3779  nbins/5,nmin,nmax/5,500,0,1.);
3780  fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
3781  fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("#it{E} (GeV)");
3782  outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
3783 
3784  if(GetCalorimeter()==kEMCAL)
3785  {
3786  fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
3787  Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
3788  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3789  fhMCEDispEta[i]->SetXTitle("#it{E} (GeV)");
3790  fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
3791  outputContainer->Add(fhMCEDispEta[i]);
3792 
3793  fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
3794  Form("cluster from %s : #sigma^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i} - <#varphi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
3795  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3796  fhMCEDispPhi[i]->SetXTitle("#it{E} (GeV)");
3797  fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3798  outputContainer->Add(fhMCEDispPhi[i]);
3799 
3800  fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
3801  Form("cluster from %s : #delta^{2}_{#eta #varphi} = #Sigma w_{i}(#varphi_{i} #eta_{i} ) / #Sigma w_{i} - <#varphi><#eta> vs E",ptypess[i].Data()),
3802  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
3803  fhMCESumEtaPhi[i]->SetXTitle("#it{E} (GeV)");
3804  fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #varphi}");
3805  outputContainer->Add(fhMCESumEtaPhi[i]);
3806 
3807  fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
3808  Form("cluster from %s : #sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),
3809  nptbins,ptmin,ptmax,200,-10,10);
3810  fhMCEDispEtaPhiDiff[i]->SetXTitle("#it{E} (GeV)");
3811  fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}-#sigma^{2}_{#eta #eta}");
3812  outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
3813 
3814  fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
3815  Form("cluster from %s : (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi}) vs E",ptypess[i].Data()),
3816  nptbins,ptmin,ptmax, 200,-1,1);
3817  fhMCESphericity[i]->SetXTitle("#it{E} (GeV)");
3818  fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi})");
3819  outputContainer->Add(fhMCESphericity[i]);
3820 
3821  for(Int_t ie = 0; ie < 7; ie++)
3822  {
3823  fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
3824  Form("cluster from %s : #sigma^{2}_{#varphi #varphi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]),
3825  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3826  fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
3827  fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3828  outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
3829 
3830  fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pnamess[i].Data()),
3831  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]),
3832  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3833  fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
3834  fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3835  outputContainer->Add(fhMCLambda0DispEta[ie][i]);
3836 
3837  fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
3838  Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#varphi #varphi} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]),
3839  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3840  fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
3841  fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3842  outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
3843  }
3844  }
3845  }
3846  }// loop
3847 
3848 // if(!GetReader()->IsEmbeddedClusterSelectionOn())
3849 // {
3850 // fhMCPhotonELambda0NoOverlap = new TH2F("hELambda0_MCPhoton_NoOverlap",
3851 // "cluster from Photon : E vs #lambda_{0}^{2}",
3852 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3853 // fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
3854 // fhMCPhotonELambda0NoOverlap->SetXTitle("#it{E} (GeV)");
3855 // outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
3856 //
3857 // fhMCPhotonELambda0TwoOverlap = new TH2F("hELambda0_MCPhoton_TwoOverlap",
3858 // "cluster from Photon : E vs #lambda_{0}^{2}",
3859 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3860 // fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
3861 // fhMCPhotonELambda0TwoOverlap->SetXTitle("#it{E} (GeV)");
3862 // outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
3863 //
3864 // fhMCPhotonELambda0NOverlap = new TH2F("hELambda0_MCPhoton_NOverlap",
3865 // "cluster from Photon : E vs #lambda_{0}^{2}",
3866 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3867 // fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
3868 // fhMCPhotonELambda0NOverlap->SetXTitle("#it{E} (GeV)");
3869 // outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
3870 // } // No embedding
3871 
3872  if(GetReader()->IsEmbeddedClusterSelectionOn())
3873  {
3874  fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
3875  "Energy Fraction of embedded signal versus cluster energy",
3876  nptbins,ptmin,ptmax,100,0.,1.);
3877  fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
3878  fhEmbeddedSignalFractionEnergy->SetXTitle("#it{E} (GeV)");
3879  outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
3880 
3881  fhEmbedPhotonELambda0FullSignal = new TH2F("hELambda0_EmbedPhoton_FullSignal",
3882  "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
3883  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3884  fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
3885  fhEmbedPhotonELambda0FullSignal->SetXTitle("#it{E} (GeV)");
3886  outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
3887 
3888  fhEmbedPhotonELambda0MostlySignal = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
3889  "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
3890  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3891  fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
3892  fhEmbedPhotonELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
3893  outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
3894 
3895  fhEmbedPhotonELambda0MostlyBkg = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
3896  "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
3897  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3898  fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
3899  fhEmbedPhotonELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
3900  outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
3901 
3902  fhEmbedPhotonELambda0FullBkg = new TH2F("hELambda0_EmbedPhoton_FullBkg",
3903  "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
3904  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3905  fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
3906  fhEmbedPhotonELambda0FullBkg->SetXTitle("#it{E} (GeV)");
3907  outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
3908 
3909  fhEmbedPi0ELambda0FullSignal = new TH2F("hELambda0_EmbedPi0_FullSignal",
3910  "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
3911  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3912  fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
3913  fhEmbedPi0ELambda0FullSignal->SetXTitle("#it{E} (GeV)");
3914  outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
3915 
3916  fhEmbedPi0ELambda0MostlySignal = new TH2F("hELambda0_EmbedPi0_MostlySignal",
3917  "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
3918  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3919  fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
3920  fhEmbedPi0ELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
3921  outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
3922 
3923  fhEmbedPi0ELambda0MostlyBkg = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
3924  "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
3925  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3926  fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
3927  fhEmbedPi0ELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
3928  outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
3929 
3930  fhEmbedPi0ELambda0FullBkg = new TH2F("hELambda0_EmbedPi0_FullBkg",
3931  "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
3932  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3933  fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
3934  fhEmbedPi0ELambda0FullBkg->SetXTitle("#it{E} (GeV)");
3935  outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
3936  }// embedded histograms
3937  }// Fill SS MC histograms
3938 
3939 
3941  {
3942  fhMCConversionVertex = new TH2F("hMCPhotonConversionVertex","cluster from converted photon, #it{p}_{T} vs vertex distance",
3943  nptbins,ptmin,ptmax,500,0,500);
3944  fhMCConversionVertex->SetYTitle("#it{R} (cm)");
3945  fhMCConversionVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3946  outputContainer->Add(fhMCConversionVertex) ;
3947 
3948  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
3949  {
3950  fhMCConversionVertexTRD = new TH2F("hMCPhotonConversionVertexTRD","cluster from converted photon, #it{p}_{T} vs vertex distance, SM covered by TRD",
3951  nptbins,ptmin,ptmax,500,0,500);
3952  fhMCConversionVertexTRD->SetYTitle("#it{R} (cm)");
3953  fhMCConversionVertexTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3954  outputContainer->Add(fhMCConversionVertexTRD) ;
3955  }
3956 
3957  if(fFillSSHistograms)
3958  {
3959  TString region[] = {"ITS","TPC","TRD","TOF","Top EMCal","In EMCal"};
3960  for(Int_t iR = 0; iR < 6; iR++)
3961  {
3962  fhMCConversionLambda0Rcut[iR] = new TH2F(Form("hMCPhotonConversionLambda0_R%d",iR),
3963  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s",region[iR].Data()),
3964  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3965  fhMCConversionLambda0Rcut[iR]->SetYTitle("#lambda_{0}^{2}");
3966  fhMCConversionLambda0Rcut[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3967  outputContainer->Add(fhMCConversionLambda0Rcut[iR]) ;
3968 
3969  fhMCConversionLambda1Rcut[iR] = new TH2F(Form("hMCPhotonConversionLambda1_R%d",iR),
3970  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{1}^{2}, conversion in %s",region[iR].Data()),
3971  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3972  fhMCConversionLambda1Rcut[iR]->SetYTitle("#lambda_{1}^{2}");
3973  fhMCConversionLambda1Rcut[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3974  outputContainer->Add(fhMCConversionLambda1Rcut[iR]) ;
3975  } // R cut
3976 
3977 
3978  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
3979  {
3980  for(Int_t iR = 0; iR < 6; iR++)
3981  {
3982  fhMCConversionLambda0RcutTRD[iR] = new TH2F(Form("hMCPhotonConversionLambda0TRD_R%d",iR),
3983  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s, SM covered by TRD",region[iR].Data()),
3984  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3985  fhMCConversionLambda0RcutTRD[iR]->SetYTitle("#lambda_{0}^{2}");
3986  fhMCConversionLambda0RcutTRD[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3987  outputContainer->Add(fhMCConversionLambda0RcutTRD[iR]) ;
3988 
3989  fhMCConversionLambda1RcutTRD[iR] = new TH2F(Form("hMCPhotonConversionLambda1TRD_R%d",iR),
3990  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{1}^{2}, conversion in %s, SM covered by TRD",region[iR].Data()),
3991  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3992  fhMCConversionLambda1RcutTRD[iR]->SetYTitle("#lambda_{1}^{2}");
3993  fhMCConversionLambda1RcutTRD[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3994  outputContainer->Add(fhMCConversionLambda1RcutTRD[iR]) ;
3995  } // R cut
3996  }
3997 
3998  // if(GetCalorimeter() == kEMCAL && fFillEMCALRegionSSHistograms)
3999  // {
4000  // for(Int_t ieta = 0; ieta < 4; ieta++)
4001  // {
4002  // for(Int_t iphi = 0; iphi < 3; iphi++)
4003  // {
4004  // for(Int_t iR = 0; iR < 6; iR++)
4005  // {
4006  // fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR] =
4007  // new TH2F(Form("hMCPhotonConversionLambda0_R%d_eta%d_phi%d",iR,ieta,iphi),
4008  // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s, region eta %d, phi %d",region[iR].Data(),ieta,iphi),
4009  // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4010  // fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR]->SetYTitle("#lambda_{0}^{2}");
4011  // fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4012  // outputContainer->Add(fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR]) ;
4013  //
4014  // if(GetFirstSMCoveredByTRD() >= 0)
4015  // {
4016  // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR] =
4017  // new TH2F(Form("hMCPhotonConversionLambda0TRD_R%d_eta%d_phi%d",iR,ieta,iphi),
4018  // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s, region eta %d, phi %d, SM covered by TRD",region[iR].Data(),ieta,iphi),
4019  // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4020  // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR]->SetYTitle("#lambda_{0}^{2}");
4021  // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4022  // outputContainer->Add(fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR]) ;
4023  // } // TRD
4024  //
4025  // } // iR
4026  // } // iphi
4027  // } // ieta
4028  // } // regions in EMCal
4029  } // shower shape
4030  } // conversion vertex
4031  } // Histos with MC
4032 
4034  {
4035  for(Int_t ie=0; ie<fNEBinCuts; ie++)
4036  {
4037  fhEBinClusterEtaPhi[ie] = new TH2F
4038  (Form("hEBin%d_Cluster_EtaPhi",ie),
4039  Form("#eta vs #varphi, cluster, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
4040  netabins,etamin,etamax,nphibins,phimin,phimax);
4041  fhEBinClusterEtaPhi[ie]->SetYTitle("#varphi (rad)");
4042  fhEBinClusterEtaPhi[ie]->SetXTitle("#eta");
4043  outputContainer->Add(fhEBinClusterEtaPhi[ie]) ;
4044 
4045  fhEBinClusterColRow[ie] = new TH2F
4046  (Form("hEBin%d_Cluster_ColRow",ie),
4047  Form("column vs row, cluster max E cell, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
4048  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
4049  fhEBinClusterColRow[ie]->SetYTitle("row");
4050  fhEBinClusterColRow[ie]->SetXTitle("column");
4051  outputContainer->Add(fhEBinClusterColRow[ie]) ;
4052 
4053  fhEBinClusterEtaPhiPID[ie] = new TH2F
4054  (Form("hEBin%d_Cluster_EtaPhi_PID",ie),
4055  Form("#eta vs #varphi, cluster, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}, PID cut",fEBinCuts[ie],fEBinCuts[ie+1]),
4056  netabins,etamin,etamax,nphibins,phimin,phimax);
4057  fhEBinClusterEtaPhiPID[ie]->SetYTitle("#varphi (rad)");
4058  fhEBinClusterEtaPhiPID[ie]->SetXTitle("#eta");
4059  outputContainer->Add(fhEBinClusterEtaPhiPID[ie]) ;
4060 
4061  fhEBinClusterColRowPID[ie] = new TH2F
4062  (Form("hEBin%d_Cluster_ColRow_PID",ie),
4063  Form("column vs row, cluster max E cell, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}, PID cut",fEBinCuts[ie],fEBinCuts[ie+1]),
4064  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
4065  fhEBinClusterColRowPID[ie]->SetYTitle("row");
4066  fhEBinClusterColRowPID[ie]->SetXTitle("column");
4067  outputContainer->Add(fhEBinClusterColRowPID[ie]) ;
4068  }
4069  }
4070 
4072  {
4073  TString caseTitle[] = {"","CleanCluster","MergedClusterHijingBkg","MergedClusterNotHijingBkg","MergedClusterHijingAndOtherBkg","MergedCluster"};
4074  Int_t ncases = 1;
4075  if ( IsDataMC() && IsStudyClusterOverlapsPerGeneratorOn() ) ncases = 6;
4076 
4077  for(Int_t icase = 0; icase < ncases; icase++)
4078  {
4079  fhLocalRegionClusterEtaPhi[icase] = new TH2F
4080  (Form("hLocalRegionClusterEtaPhi%s",caseTitle[icase].Data()),
4081  "cluster,#it{E} > 0.5 GeV, #eta vs #varphi",netabins,etamin,etamax,nphibins,phimin,phimax);
4082  fhLocalRegionClusterEtaPhi[icase]->SetYTitle("#varphi (rad)");
4083  fhLocalRegionClusterEtaPhi[icase]->SetXTitle("#eta");
4084  outputContainer->Add(fhLocalRegionClusterEtaPhi[icase]) ;
4085 
4086  fhLocalRegionClusterEnergySum[icase] = new TH2F
4087  (Form("hLocalRegionClusterEnergySum%s",caseTitle[icase].Data()),
4088  "Sum of cluster energy around trigger cluster #it{E} with R=0.2",
4089  nptbins,ptmin,ptmax, 200,0,100);
4090  fhLocalRegionClusterEnergySum[icase]->SetXTitle("#it{E} (GeV)");
4091  fhLocalRegionClusterEnergySum[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4092  outputContainer->Add(fhLocalRegionClusterEnergySum[icase]);
4093 
4095  (Form("hLocalRegionClusterMultiplicity%s",caseTitle[icase].Data()),
4096  "Cluster multiplicity around trigger cluster #it{E} with R=0.2",
4097  nptbins,ptmin,ptmax, 200,0,200);
4098  fhLocalRegionClusterMultiplicity[icase]->SetXTitle("#it{E} (GeV)");
4099  fhLocalRegionClusterMultiplicity[icase]->SetYTitle("Multiplicity");
4100  outputContainer->Add(fhLocalRegionClusterMultiplicity[icase]);
4101 
4103  {
4105  (Form("hLocalRegionClusterEnergySumPerCentrality%s",caseTitle[icase].Data()),
4106  "Sum of cluster energy around trigger cluster vs centrality with R=0.2",
4107  100,0,100, 200,0,100);
4108  fhLocalRegionClusterEnergySumPerCentrality[icase]->SetXTitle("Centrality");
4109  fhLocalRegionClusterEnergySumPerCentrality[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4110  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentrality[icase]);
4111 
4113  (Form("hLocalRegionClusterMultiplicityPerCentrality%s",caseTitle[icase].Data()),
4114  "Cluster multiplicity around trigger cluster vs centrality with R=0.2",
4115  100,0,100, 200,0,200);
4116  fhLocalRegionClusterMultiplicityPerCentrality[icase]->SetXTitle("Centrality");
4117  fhLocalRegionClusterMultiplicityPerCentrality[icase]->SetYTitle("Multiplicity");
4118  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentrality[icase]);
4119  }
4120 
4121  //
4122  // Select only single pi0 decay clusters from non hijing generator
4123  //
4124  if(IsDataMC())
4125  {
4127  (Form("hLocalRegionClusterEnergySum%s_MCPi0Decay",caseTitle[icase].Data()),
4128  "Sum of cluster energy around trigger cluster #it{E} with R=0.2",
4129  nptbins,ptmin,ptmax, 200,0,100);
4130  fhLocalRegionClusterEnergySumMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4131  fhLocalRegionClusterEnergySumMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4132  outputContainer->Add(fhLocalRegionClusterEnergySumMCPi0Decay[icase]);
4133 
4135  (Form("hLocalRegionClusterMultiplicity%s_MCPi0Decay",caseTitle[icase].Data()),
4136  "Cluster multiplicity around trigger cluster #it{E} with R=0.2",
4137  nptbins,ptmin,ptmax, 200,0,200);
4138  fhLocalRegionClusterMultiplicityMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4139  fhLocalRegionClusterMultiplicityMCPi0Decay[icase]->SetYTitle("Multiplicity");
4140  outputContainer->Add(fhLocalRegionClusterMultiplicityMCPi0Decay[icase]);
4141 
4143  {
4145  (Form("hLocalRegionClusterEnergySumPerCentrality%s_MCPi0Decay",caseTitle[icase].Data()),
4146  "Sum of cluster energy around trigger cluster vs centrality with R=0.2",
4147  100,0,100, 200,0,100);
4148  fhLocalRegionClusterEnergySumPerCentralityMCPi0Decay[icase]->SetXTitle("Centrality");
4149  fhLocalRegionClusterEnergySumPerCentralityMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4150  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityMCPi0Decay[icase]);
4151 
4153  (Form("hLocalRegionClusterMultiplicityPerCentrality%s_MCPi0Decay",caseTitle[icase].Data()),
4154  "Cluster multiplicity around trigger cluster vs centrality with R=0.2",
4155  100,0,100, 200,0,200);
4156  fhLocalRegionClusterMultiplicityPerCentralityMCPi0Decay[icase]->SetXTitle("Centrality");
4157  fhLocalRegionClusterMultiplicityPerCentralityMCPi0Decay[icase]->SetYTitle("Multiplicity");
4158  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentralityMCPi0Decay[icase]);
4159  }
4160  }
4161 
4163  {
4165  (Form("hLocalRegionClusterEnergySumHijing%s",caseTitle[icase].Data()),
4166  "Sum of cluster energy (HIJING) around trigger cluster #it{E} with R=0.2",
4167  nptbins,ptmin,ptmax, 200,0,100);
4168  fhLocalRegionClusterEnergySumHijing[icase]->SetXTitle("#it{E} (GeV)");
4169  fhLocalRegionClusterEnergySumHijing[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4170  outputContainer->Add(fhLocalRegionClusterEnergySumHijing[icase]);
4171 
4173  (Form("hLocalRegionClusterMultiplicityHijing%s",caseTitle[icase].Data()),
4174  "Cluster multiplicity (HIJING) around trigger cluster #it{E} with R=0.2",
4175  nptbins,ptmin,ptmax, 200,0,200);
4176  fhLocalRegionClusterMultiplicityHijing[icase]->SetXTitle("#it{E} (GeV)");
4177  fhLocalRegionClusterMultiplicityHijing[icase]->SetYTitle("Multiplicity");
4178  outputContainer->Add(fhLocalRegionClusterMultiplicityHijing[icase]);
4179 
4181  (Form("hLocalRegionClusterEnergySumAdded%s",caseTitle[icase].Data()),
4182  "Sum of cluster energy (not HIJING) around trigger cluster #it{E} with R=0.2",
4183  nptbins,ptmin,ptmax, 200,0,100);
4184  fhLocalRegionClusterEnergySumAdded[icase]->SetXTitle("#it{E} (GeV)");
4185  fhLocalRegionClusterEnergySumAdded[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4186  outputContainer->Add(fhLocalRegionClusterEnergySumAdded[icase]);
4187 
4189  (Form("hLocalRegionClusterMultiplicityAdded%s",caseTitle[icase].Data()),
4190  "Cluster multiplicity (not HIJING) around trigger cluster #it{E} with R=0.2",
4191  nptbins,ptmin,ptmax, 200,0,200);
4192  fhLocalRegionClusterMultiplicityAdded[icase]->SetXTitle("#it{E} (GeV)");
4193  fhLocalRegionClusterMultiplicityAdded[icase]->SetYTitle("Multiplicity");
4194  outputContainer->Add(fhLocalRegionClusterMultiplicityAdded[icase]);
4195 
4196 
4198  {
4200  (Form("hLocalRegionClusterEnergySumPerCentralityHijing%s",caseTitle[icase].Data()),
4201  "Sum of cluster energy (HIJING) around trigger cluster vs centrality with R=0.2",
4202  100,0,100, 200,0,100);
4203  fhLocalRegionClusterEnergySumPerCentralityHijing[icase]->SetXTitle("Centrality");
4204  fhLocalRegionClusterEnergySumPerCentralityHijing[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4205  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityHijing[icase]);
4206 
4208  (Form("hLocalRegionClusterMultiplicityPerCentralityHijing%s",caseTitle[icase].Data()),
4209  "Cluster multiplicity (HIJING) around trigger cluster vs centrality with R=0.2",
4210  100,0,100, 200,0,200);
4211  fhLocalRegionClusterMultiplicityPerCentralityHijing[icase]->SetXTitle("Centrality");
4212  fhLocalRegionClusterMultiplicityPerCentralityHijing[icase]->SetYTitle("Multiplicity");
4213  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentralityHijing[icase]);
4214 
4215 
4217  (Form("hLocalRegionClusterEnergySumPerCentralityAdded%s",caseTitle[icase].Data()),
4218  "Sum of cluster energy (not HIJING) around trigger cluster vs centrality with R=0.2",
4219  100,0,100, 200,0,100);
4220  fhLocalRegionClusterEnergySumPerCentralityAdded[icase]->SetXTitle("Centrality");
4221  fhLocalRegionClusterEnergySumPerCentralityAdded[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4222  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityAdded[icase]);
4223 
4225  (Form("hLocalRegionClusterMultiplicityPerCentralityAdded%s",caseTitle[icase].Data()),
4226  "Cluster multiplicity (not HIJING) around trigger cluster vs centrality with R=0.2",
4227  100,0,100, 200,0,200);
4228  fhLocalRegionClusterMultiplicityPerCentralityAdded[icase]->SetXTitle("Centrality");
4229  fhLocalRegionClusterMultiplicityPerCentralityAdded[icase]->SetYTitle("Multiplicity");
4230  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentralityAdded[icase]);
4231  }
4232 
4233  //
4234  // Select only single pi0 decay clusters from non hijing generator
4235  //
4236 
4238  (Form("hLocalRegionClusterEnergySumHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4239  "Sum of cluster energy (HIJING) around trigger cluster #it{E} with R=0.2",
4240  nptbins,ptmin,ptmax, 200,0,100);
4241  fhLocalRegionClusterEnergySumHijingMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4242  fhLocalRegionClusterEnergySumHijingMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4243  outputContainer->Add(fhLocalRegionClusterEnergySumHijingMCPi0Decay[icase]);
4244 
4246  (Form("hLocalRegionClusterMultiplicityHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4247  "Cluster multiplicity (HIJING) around trigger cluster #it{E} with R=0.2",
4248  nptbins,ptmin,ptmax, 200,0,200);
4249  fhLocalRegionClusterMultiplicityHijingMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4250  fhLocalRegionClusterMultiplicityHijingMCPi0Decay[icase]->SetYTitle("Multiplicity");
4251  outputContainer->Add(fhLocalRegionClusterMultiplicityHijingMCPi0Decay[icase]);
4252 
4254  (Form("hLocalRegionClusterEnergySumAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4255  "Sum of cluster energy (not HIJING) around trigger cluster #it{E} with R=0.2",
4256  nptbins,ptmin,ptmax, 200,0,100);
4257  fhLocalRegionClusterEnergySumAddedMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4258  fhLocalRegionClusterEnergySumAddedMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4259  outputContainer->Add(fhLocalRegionClusterEnergySumAddedMCPi0Decay[icase]);
4260 
4262  (Form("hLocalRegionClusterMultiplicityAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4263  "Cluster multiplicity (not HIJING) around trigger cluster #it{E} with R=0.2",
4264  nptbins,ptmin,ptmax, 200,0,200);
4265  fhLocalRegionClusterMultiplicityAddedMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4266  fhLocalRegionClusterMultiplicityAddedMCPi0Decay[icase]->SetYTitle("Multiplicity");
4267  outputContainer->Add(fhLocalRegionClusterMultiplicityAddedMCPi0Decay[icase]);
4268 
4269 
4271  {
4273  (Form("hLocalRegionClusterEnergySumPerCentralityHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4274  "Sum of cluster energy (HIJING) around trigger cluster vs centrality with R=0.2",
4275  100,0,100, 200,0,100);
4276  fhLocalRegionClusterEnergySumPerCentralityHijingMCPi0Decay[icase]->SetXTitle("Centrality");
4277  fhLocalRegionClusterEnergySumPerCentralityHijingMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4278  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityHijingMCPi0Decay[icase]);
4279 
4281  (Form("hLocalRegionClusterMultiplicityPerCentralityHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4282  "Cluster multiplicity (HIJING) around trigger cluster vs centrality with R=0.2",
4283  100,0,100, 200,0,200);
4284  fhLocalRegionClusterMultiplicityPerCentralityHijingMCPi0Decay[icase]->SetXTitle("Centrality");
4285  fhLocalRegionClusterMultiplicityPerCentralityHijingMCPi0Decay[icase]->SetYTitle("Multiplicity");
4287 
4288 
4290  (Form("hLocalRegionClusterEnergySumPerCentralityAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4291  "Sum of cluster energy (not HIJING) around trigger cluster vs centrality with R=0.2",
4292  100,0,100, 200,0,100);
4293  fhLocalRegionClusterEnergySumPerCentralityAddedMCPi0Decay[icase]->SetXTitle("Centrality");
4294  fhLocalRegionClusterEnergySumPerCentralityAddedMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4295  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityAddedMCPi0Decay[icase]);
4296 
4298  (Form("hLocalRegionClusterMultiplicityPerCentralityAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4299  "Cluster multiplicity (not HIJING) around trigger cluster vs centrality with R=0.2",
4300  100,0,100, 200,0,200);
4301  fhLocalRegionClusterMultiplicityPerCentralityAddedMCPi0Decay[icase]->SetXTitle("Centrality");
4302  fhLocalRegionClusterMultiplicityPerCentralityAddedMCPi0Decay[icase]->SetYTitle("Multiplicity");
4304  }
4305 
4306  }
4307  }
4308 
4309  if(IsDataMC())
4310  {
4312  ("hLocalRegionClusterEnergySumHijing2",
4313  "Sum of cluster energy (HIJING) around trigger cluster #it{E} with R=0.2",
4314  nptbins,ptmin,ptmax, 200,0,100);
4315  fhLocalRegionClusterEnergySumHijing2->SetXTitle("#it{E} (GeV)");
4316  fhLocalRegionClusterEnergySumHijing2->SetYTitle("#Sigma #it{E} (GeV)");
4317  outputContainer->Add(fhLocalRegionClusterEnergySumHijing2);
4318 
4320  ("hLocalRegionClusterMultiplicityHijing2",
4321  "Cluster multiplicity (HIJING) around trigger cluster #it{E} with R=0.2",
4322  nptbins,ptmin,ptmax, 200,0,200);
4323  fhLocalRegionClusterMultiplicityHijing2->SetXTitle("#it{E} (GeV)");
4324  fhLocalRegionClusterMultiplicityHijing2->SetYTitle("Multiplicity");
4325  outputContainer->Add(fhLocalRegionClusterMultiplicityHijing2);
4326 
4328  {
4330  ("hLocalRegionClusterEnergySumPerCentralityHijing2",
4331  "Sum of cluster energy (HIJING) around trigger cluster vs centrality with R=0.2",
4332  100,0,100, 200,0,100);
4333  fhLocalRegionClusterEnergySumPerCentralityHijing2->SetXTitle("Centrality");
4334  fhLocalRegionClusterEnergySumPerCentralityHijing2->SetYTitle("#Sigma #it{E} (GeV)");
4336 
4338  ("hLocalRegionClusterMultiplicityPerCentralityHijing2",
4339  "Cluster multiplicity (HIJING) around trigger cluster vs centrality with R=0.2",
4340  100,0,100, 200,0,200);
4342  fhLocalRegionClusterMultiplicityPerCentralityHijing2->SetYTitle("Multiplicity");
4344  }
4345  }
4346  // fhDistanceAddedPhotonAddedPrimarySignal = new TH2F
4347  // ("hDistanceAddedPhotonAddedPrimarySignal", "Distance added #gamma to primary particle from added generator"
4348  // ,nptbins,ptmin,ptmax,100,0,0.4);
4349  // fhDistanceAddedPhotonAddedPrimarySignal->SetYTitle("#it{R}");
4350  // fhDistanceAddedPhotonAddedPrimarySignal->SetXTitle("#it{p}_{T}");
4351  // outputContainer->Add(fhDistanceAddedPhotonAddedPrimarySignal) ;
4352  //
4353  // fhDistanceHijingPhotonAddedPrimarySignal = new TH2F
4354  // ("hDistanceHijingPhotonAddedPrimarySignal", "Distance Hijing #gamma to primary particle from added generator"
4355  // ,nptbins,ptmin,ptmax,100,0,0.4);
4356  // fhDistanceHijingPhotonAddedPrimarySignal->SetYTitle("#it{R}");
4357  // fhDistanceHijingPhotonAddedPrimarySignal->SetXTitle("#it{p}_{T}");
4358  // outputContainer->Add(fhDistanceHijingPhotonAddedPrimarySignal) ;
4359  //
4360  // fhDistanceAddedPhotonAddedSecondarySignal = new TH2F
4361  // ("hDistanceAddedPhotonAddedSecondarySignal", "Distance added #gamma to secondary particle from added generator"
4362  // ,nptbins,ptmin,ptmax,100,0,0.4);
4363  // fhDistanceAddedPhotonAddedSecondarySignal->SetYTitle("#it{R}");
4364  // fhDistanceAddedPhotonAddedSecondarySignal->SetXTitle("#it{p}_{T}");
4365  // outputContainer->Add(fhDistanceAddedPhotonAddedSecondarySignal) ;
4366  //
4367  // fhDistanceHijingPhotonAddedSecondarySignal = new TH2F
4368  // ("hDistanceHijingPhotonAddedSecondarySignal", "Distance Hijing #gamma to secondary particle from added generator"
4369  // ,nptbins,ptmin,ptmax,100,0,0.4);
4370  // fhDistanceHijingPhotonAddedSecondarySignal->SetYTitle("#it{R}");
4371  // fhDistanceHijingPhotonAddedSecondarySignal->SetXTitle("#it{p}_{T}");
4372  // outputContainer->Add(fhDistanceHijingPhotonAddedSecondarySignal) ;
4373  //
4374  // fhDistanceAddedPhotonHijingSecondary = new TH2F
4375  // ("hDistanceAddedPhotonHijingSecondary", "Distance added #gamma to secondary particle from hijing"
4376  // ,nptbins,ptmin,ptmax,100,0,0.4);
4377  // fhDistanceAddedPhotonHijingSecondary->SetYTitle("#it{R}");
4378  // fhDistanceAddedPhotonHijingSecondary->SetXTitle("#it{p}_{T}");
4379  // outputContainer->Add(fhDistanceAddedPhotonHijingSecondary) ;
4380 
4381 
4383  ("hDistance2AddedSignals", "Distance added signals"
4384  ,nptbins,ptmin,ptmax,100,0,0.4);
4385  fhDistance2AddedSignals->SetYTitle("#it{R}");
4386  fhDistance2AddedSignals->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4387  outputContainer->Add(fhDistance2AddedSignals) ;
4388 
4389  fhDistance2Hijing = new TH2F
4390  ("hDistance2Hijing", "Distance 2 hijing clusters"
4391  ,nptbins,ptmin,ptmax,100,0,0.4);
4392  fhDistance2Hijing->SetYTitle("#it{R}");
4393  fhDistance2Hijing->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4394  outputContainer->Add(fhDistance2Hijing) ;
4395 
4397  ("hDistanceAddedSignalsHijing", "Distance added signals to hijing"
4398  ,nptbins,ptmin,ptmax,100,0,0.4);
4399  fhDistanceAddedSignalsHijing->SetYTitle("#it{R}");
4400  fhDistanceAddedSignalsHijing->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4401  outputContainer->Add(fhDistanceAddedSignalsHijing) ;
4402  }
4403 
4405  {
4406 
4407  TString mcGenNames[] = {"","_MC_Pi0Merged","_MC_Pi0Decay","_MC_EtaDecay","_MC_PhotonOther","_MC_Electron","_MC_Other"};
4408  TString mcGenTitle[] = {"",",MC Pi0-Merged",",MC Pi0-Decay",", MC Eta-Decay",", MC Photon other sources",", MC Electron",", MC other sources"};
4409  for(Int_t igen = 0; igen < GetNCocktailGenNamesToCheck(); igen++)
4410  {
4411  TString add = "_MainGener_";
4412  if(igen==0) add = "";
4413  for(Int_t imc = 0; imc < fgkNGenTypes; imc++)
4414  {
4415  fhCleanGeneratorCluster[igen][imc] = new TH1F(Form("hCleanGeneratorCluster%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4416  Form("Number of selected clusters with contribution of %s generator%s, no overlap",
4417  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4418  nptbins,ptmin,ptmax);
4419  fhCleanGeneratorCluster[igen][imc]->SetYTitle("#it{counts}");
4420  fhCleanGeneratorCluster[igen][imc]->SetXTitle("#it{E} (GeV)");
4421  outputContainer->Add(fhCleanGeneratorCluster[igen][imc]) ;
4422 
4423  fhCleanGeneratorClusterEPrimRecoRatio[igen][imc] = new TH2F(Form("hCleanGeneratorClusterEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4424  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of %s generator%s, no overlap",
4425  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4426  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4427  fhCleanGeneratorClusterEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4428  fhCleanGeneratorClusterEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4429  outputContainer->Add(fhCleanGeneratorClusterEPrimRecoRatio[igen][imc]) ;
4430 
4431  fhCleanGeneratorClusterEPrimRecoDiff[igen][imc] = new TH2F(Form("hCleanGeneratorClusterEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4432  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of %s generator%s, no overlap",
4433  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4434  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4435  fhCleanGeneratorClusterEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4436  fhCleanGeneratorClusterEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4437  outputContainer->Add(fhCleanGeneratorClusterEPrimRecoDiff[igen][imc]) ;
4438 
4439  //
4440 
4441  fhMergeGeneratorCluster[igen][imc] = new TH1F(Form("hMergeGeneratorCluster%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4442  Form("Number of selected clusters with contribution of >=2 generators, main %s%s",
4443  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4444  nptbins,ptmin,ptmax);
4445  fhMergeGeneratorCluster[igen][imc]->SetYTitle("#it{counts}");
4446  fhMergeGeneratorCluster[igen][imc]->SetXTitle("#it{E} (GeV)");
4447  outputContainer->Add(fhMergeGeneratorCluster[igen][imc]) ;
4448 
4449  fhMergeGeneratorClusterEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4450  Form("#it{E}_{reco}/#it{E}_{gen}clusters with contribution of >=2 generators, main %s%s",
4451  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4452  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4453  fhMergeGeneratorClusterEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4454  fhMergeGeneratorClusterEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4455  outputContainer->Add(fhMergeGeneratorClusterEPrimRecoRatio[igen][imc]) ;
4456 
4457  fhMergeGeneratorClusterEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4458  Form("#it{E}_{reco}-#it{E}_{gen}clusters with contribution of >=2 generators, main %s%s",
4459  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4460  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4461  fhMergeGeneratorClusterEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4462  fhMergeGeneratorClusterEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4463  outputContainer->Add(fhMergeGeneratorClusterEPrimRecoDiff[igen][imc]) ;
4464 
4465  //if(GetCocktailGenNameToCheck(igen).Contains("ijing")) continue;
4466 
4467  //
4468 
4469  fhMergeGeneratorClusterNotHijingBkg[igen][imc] = new TH1F(Form("hMergeGeneratorClusterNotHijingBkg%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4470  Form("Number of selected clusters with contribution of >=2 generators, , none is HIJING, main %s%s",
4471  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4472  nptbins,ptmin,ptmax);
4473  fhMergeGeneratorClusterNotHijingBkg[igen][imc]->SetYTitle("#it{counts}");
4474  fhMergeGeneratorClusterNotHijingBkg[igen][imc]->SetXTitle("#it{E} (GeV)");
4475  outputContainer->Add(fhMergeGeneratorClusterNotHijingBkg[igen][imc]) ;
4476 
4477  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterNotHijingBkgEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4478  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of >=2 generators, , none is HIJING, main %s%s",
4479  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4480  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4481  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4482  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4483  outputContainer->Add(fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc]) ;
4484 
4485  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterNotHijingBkgEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4486  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of >=2 generators, , none is HIJING, main %s%s",
4487  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4488  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4489  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4490  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4491  outputContainer->Add(fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc]) ;
4492 
4493  //
4494 
4495  fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc] = new TH1F(Form("hMergeGeneratorClusterHijingAndOtherBkg%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4496  Form("Number of selected clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4497  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4498  nptbins,ptmin,ptmax);
4499  fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc]->SetYTitle("#it{counts}");
4500  fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc]->SetXTitle("#it{E} (GeV)");
4501  outputContainer->Add(fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc]) ;
4502 
4503 
4504  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4505  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4506  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4507  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4508  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4509  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4510  outputContainer->Add(fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc]) ;
4511 
4512 
4513  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4514  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4515  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4516  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4517  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4518  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4519  outputContainer->Add(fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc]) ;
4520 
4521  //
4522 
4523  fhMergeGeneratorClusterHijingBkg[igen][imc] = new TH1F(Form("hMergeGeneratorClusterHijingBkg%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4524  Form("Number of selected clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4525  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4526  nptbins,ptmin,ptmax);
4527  fhMergeGeneratorClusterHijingBkg[igen][imc]->SetYTitle("#it{counts}");
4528  fhMergeGeneratorClusterHijingBkg[igen][imc]->SetXTitle("#it{E} (GeV)");
4529  outputContainer->Add(fhMergeGeneratorClusterHijingBkg[igen][imc]) ;
4530 
4531  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingBkgEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4532  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4533  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4534  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4535  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4536  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4537  outputContainer->Add(fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc]) ;
4538 
4539  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingBkgEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4540  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4541  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4542  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4543  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4544  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4545  outputContainer->Add(fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc]) ;
4546  }
4547  }
4548  }
4549 
4550  return outputContainer ;
4551 }
4552 
4553 //_______________________
4555 //_______________________
4557 {
4558  if ( GetCalorimeter() == kPHOS && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD() )
4559  AliFatal("!!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
4560  else if( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
4561  AliFatal("!!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
4562 
4563  // Trick to select primary photon particles in the analysis of pure MC input
4564  if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
4565 }
4566 
4567 //_________________________________
4569 //_________________________________
4571 {
4572  AddToHistogramsName("AnaPhoton_");
4573 
4574  fMinDist = 2.;
4575  fMinDist2 = 4.;
4576  fMinDist3 = 5.;
4577 
4578  fTimeCutMin =-1000000;
4579  fTimeCutMax = 1000000;
4580  fNCellsCut = 0;
4581 
4582  fRejectTrackMatch = kTRUE ;
4583 
4584  fNEBinCuts = 14;
4585  fEBinCuts[0] = 0.; fEBinCuts[1] = 0.3; fEBinCuts[2] = 0.5;
4586  fEBinCuts[3] = 1.; fEBinCuts[4] = 2. ; fEBinCuts[5] = 3. ;
4587  fEBinCuts[6] = 4.; fEBinCuts[7] = 5. ; fEBinCuts[8] = 7. ;
4588  fEBinCuts[9] = 9.; fEBinCuts[10]= 12.; fEBinCuts[11]= 15.;
4589  fEBinCuts[12]= 20.; fEBinCuts[13]= 50.; fEBinCuts[14]= 100.;
4590  for(Int_t i = fNEBinCuts; i < 15; i++) fEBinCuts[i] = 1000.;
4591 }
4592 
4593 //_______________________________________
4597 //_______________________________________
4599 {
4600  // Get the vertex
4601  Double_t v[3] = {0,0,0}; //vertex ;
4602  GetReader()->GetVertex(v);
4603 
4604  // Select the Calorimeter of the photon
4605  TObjArray * pl = 0x0;
4606  AliVCaloCells* cells = 0;
4607  if (GetCalorimeter() == kPHOS )
4608  {
4609  pl = GetPHOSClusters();
4610  cells = GetPHOSCells();
4611  }
4612  else if (GetCalorimeter() == kEMCAL)
4613  {
4614  pl = GetEMCALClusters();
4615  cells = GetEMCALCells();
4616  }
4617 
4618  if(!pl)
4619  {
4620  AliWarning(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
4621  return;
4622  }
4623 
4624  // Loop on raw clusters before filtering in the reader and fill control histogram
4625  if((GetReader()->GetEMCALClusterListName()=="" && GetCalorimeter()==kEMCAL) || GetCalorimeter()==kPHOS)
4626  {
4627  for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
4628  {
4629  AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
4630 
4631  if (GetCalorimeter() == kPHOS && clus->IsPHOS() && clus->E() > GetReader()->GetPHOSPtMin() )
4632  {
4633  clus->GetMomentum(fMomentum,GetVertex(0)) ;
4634 
4635  fhClusterCutsE [0]->Fill(fMomentum.E() , GetEventWeight());
4636  fhClusterCutsPt[0]->Fill(fMomentum.Pt(), GetEventWeight());
4637  }
4638  else if(GetCalorimeter() == kEMCAL && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin())
4639  {
4640  clus->GetMomentum(fMomentum,GetVertex(0)) ;
4641 
4642  fhClusterCutsE [0]->Fill(fMomentum.E(), GetEventWeight());
4643  fhClusterCutsPt[0]->Fill(fMomentum.Pt(), GetEventWeight());
4644  }
4645  }
4646  }
4647  else
4648  { // reclusterized
4649  TClonesArray * clusterList = 0;
4650 
4651  if(GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()))
4652  clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
4653  else if(GetReader()->GetOutputEvent())
4654  clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
4655 
4656  if(clusterList)
4657  {
4658  Int_t nclusters = clusterList->GetEntriesFast();
4659  for (Int_t iclus = 0; iclus < nclusters; iclus++)
4660  {
4661  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
4662 
4663  if(clus && clus->E() > GetReader()->GetEMCALPtMin())
4664  {
4665  clus->GetMomentum(fMomentum,GetVertex(0)) ;
4666 
4667  fhClusterCutsE [0]->Fill(clus->E() , GetEventWeight());
4668  fhClusterCutsPt[0]->Fill(fMomentum.Pt(), GetEventWeight());
4669  }
4670  }
4671  }
4672  }
4673 
4674  // Init arrays, variables, get number of clusters
4675  Int_t nCaloClusters = pl->GetEntriesFast();
4676 
4677  AliDebug(1,Form("Input %s cluster entries %d", GetCalorimeterString().Data(), nCaloClusters));
4678 
4679  //----------------------------------------------------
4680  // Fill AOD with PHOS/EMCAL AliCaloTrackParticle objects
4681  //----------------------------------------------------
4682  // Loop on clusters
4683  for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
4684  {
4685  AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
4686  //printf("calo %d, %f\n",icalo,calo->E());
4687 
4688  //Get the index where the cluster comes, to retrieve the corresponding vertex
4689  Int_t evtIndex = 0 ;
4690  if (GetMixedEvent())
4691  {
4692  evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
4693  //Get the vertex and check it is not too large in z
4694  if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
4695  }
4696 
4697  //Cluster selection, not charged, with photon id and in fiducial cut
4699  {
4700  calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
4701  }//Assume that come from vertex in straight line
4702  else
4703  {
4704  Double_t vertex[]={0,0,0};
4705  calo->GetMomentum(fMomentum,vertex) ;
4706  }
4707 
4708  //-------------------------------------
4709  // Play with the MC stack if available
4710  //-------------------------------------
4711 
4712  // Check origin of the candidates
4713  // recover different MC information
4714  Bool_t ok = kFALSE;
4715  Int_t pdg = 0, status = 0, momLabel = -1;
4716  Int_t tag = -1;
4717  Float_t mcbin = -1;
4718  Int_t noverlaps = -1;
4719  Float_t egen = 0;
4720  Bool_t conversion = kFALSE;
4721 
4722  const Int_t mcLabel = calo->GetLabel();
4723  const Int_t nlabels = calo->GetNLabels();
4724  Float_t ener = calo->E();
4725 
4726  if ( IsDataMC() && mcLabel >= 0 )
4727  {
4728  tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),
4729  calo->GetClusterMCEdepFraction(),
4730  nlabels,
4731  GetMC(),
4733  fMomentum.E(),
4734  pl); // check lost decays
4735 
4736  AliDebug(1,Form("Origin of candidate, bit map %d",tag));
4737 
4739 
4740 // GetMCAnalysisUtils()->PrintMCTag(tag);
4741 // GetMCAnalysisUtils()->PrintAncestry(GetMC(),calo->GetLabel());
4742 
4743  fPrimaryMom = GetMCAnalysisUtils()->GetMother(mcLabel, GetMC(), pdg, status, ok, momLabel);
4744 
4745  // In some histograms, each particle is a different bin
4746  if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ){ mcbin = 1.5; pdg = 111;
4747  fPrimaryMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel, 111, GetMC(),ok, momLabel); }
4748  else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ){ mcbin = 1.5; pdg = 221;
4749