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