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