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