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