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