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