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