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