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