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