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