AliPhysics  cbd2963 (cbd2963)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaPhoton.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TH2F.h>
18 #include <TClonesArray.h>
19 #include <TObjString.h>
20 #include "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 
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), 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("ESD 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("ESD 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 
1391  Float_t pt = fMomentum.Pt();
1392  Float_t eta = fMomentum.Eta();
1393  Float_t phi = fMomentum.Phi();
1394  if(phi < 0) phi+=TMath::TwoPi();
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  fhDispE ->Fill(energy, disp , GetEventWeight());
1404  fhDispPt->Fill(pt , disp , GetEventWeight());
1405  }
1406 
1408  {
1409  if(nlm==1)
1410  {
1411  fhLam0PtNLM1->Fill(pt, lambda0, GetEventWeight());
1412  fhLam1PtNLM1->Fill(pt, lambda1, GetEventWeight());
1413  }
1414  else if(nlm==2)
1415  {
1416  fhLam0PtNLM2->Fill(pt, lambda0, GetEventWeight());
1417  fhLam1PtNLM2->Fill(pt, lambda1, GetEventWeight());
1418  }
1419  }
1420 
1421  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
1422  GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
1423  {
1424  fhLam0ETRD ->Fill(energy, lambda0, GetEventWeight());
1425  fhLam0PtTRD->Fill(pt , lambda0, GetEventWeight());
1426  fhLam1ETRD ->Fill(energy, lambda1, GetEventWeight());
1428  fhDispETRD ->Fill(energy, disp, GetEventWeight());
1429  }
1430 
1431  //
1432  // EMCal SM regions
1433  //
1434  if(cluster->IsEMCAL() && fFillEMCALRegionSSHistograms)
1435  {
1436  Int_t sm = GetModuleNumber(cluster);
1437 
1438 // Bool_t shared = GetCaloUtils()->IsClusterSharedByTwoSuperModules(GetEMCALGeometry(),cluster);
1439 
1440  fhLam0PerSM[sm]->Fill(pt, lambda0, GetEventWeight());
1441  fhLam1PerSM[sm]->Fill(pt, lambda1, GetEventWeight());
1442 
1443 // if(GetReader()->IsPileUpFromSPD())
1444 // {
1445 // fhLam0PerSMSPDPileUp[sm]->Fill(pt, lambda0, GetEventWeight());
1446 // fhLam1PerSMSPDPileUp[sm]->Fill(pt, lambda1, GetEventWeight());
1447 // }
1448 
1449  Int_t etaRegion = -1, phiRegion = -1;
1450  GetCaloUtils()->GetEMCALSubregion(cluster,GetReader()->GetEMCALCells(),etaRegion,phiRegion);
1451  if(etaRegion >= 0 && etaRegion < 4 && phiRegion >=0 && phiRegion < 3)
1452  {
1453 
1454  // fhLam0EMCALRegion[etaRegion][phiRegion]->Fill(pt,lambda0, GetEventWeight());
1455  //
1456  // if(GetFirstSMCoveredByTRD() >= 0 && sm >= GetFirstSMCoveredByTRD() )
1457  // {
1458  // fhLam0EMCALRegionTRD[etaRegion][phiRegion]->Fill(pt, lambda0, GetEventWeight());
1459  // }
1460 
1461  fhLam0EMCALRegionPerSM[etaRegion][phiRegion][sm]->Fill(pt,lambda0, GetEventWeight());
1462  fhLam1EMCALRegionPerSM[etaRegion][phiRegion][sm]->Fill(pt,lambda1, GetEventWeight());
1463 
1464 // if(shared)
1465 // {
1466 // fhLam0PerSMShared[sm]->Fill(pt,lambda0, GetEventWeight());
1467 // fhLam1PerSMShared[sm]->Fill(pt,lambda1, GetEventWeight());
1468 // }
1469  }
1470 
1471  // Define l0 bin and cluster pt bin
1472  Int_t l0bin = -1;
1473  if ( lambda0 >=0.30 && lambda0 <= 0.40 ) l0bin = 1;
1474  else if( lambda0 >=0.23 && lambda0 <= 0.26 ) l0bin = 0;
1475 
1476  Float_t ptLimit[] = {2,3,4,5,6,8,10,12};
1477  Int_t ptbin = -1;
1478  for(Int_t ipt = 0; ipt < 7; ipt++)
1479  {
1480  if( pt >= ptLimit[ipt] && pt < ptLimit[ipt+1] )
1481  {
1482  ptbin = ipt;
1483  break;
1484  }
1485  }
1486 
1487  if(l0bin!=-1)
1488  {
1489  fhLam1Lam0BinPerSM[l0bin][sm]->Fill(pt,lambda1);
1490 
1491  if(ptbin >= 0)
1492  {
1493  fhEtaPhiLam0BinPtBin[l0bin][ptbin]->Fill(eta, phi, GetEventWeight());
1494 // if(shared) fhEtaPhiLam0BinPtBinSMShared[l0bin][ptbin]->Fill(eta, phi, GetEventWeight());
1495  }
1496  }
1497 
1498  //
1499  // Histograms with time of secondary cells, only if they contribute with a non null
1500  // weight to the shower shape
1501  //
1502  AliVCaloCells* cells = GetReader()->GetEMCALCells();
1503 
1504  Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
1505 
1506  //
1507  // Sort the cluster cells in energy, find maximum energy cell
1508  //
1509  Int_t ncell = cluster->GetNCells();
1510  Int_t * sortList = new Int_t [ncell];
1511  Float_t * enerList = new Float_t [ncell];
1512  Double_t * timeList = new Double_t[ncell];
1513 
1514  for(Int_t icell = 0; icell < ncell; icell++)
1515  {
1516  Int_t absId = cluster->GetCellAbsId(icell);
1517  Float_t cellE = cells->GetCellAmplitude(absId);
1518  Double_t cellTime = cells->GetCellTime(absId);
1519 
1520  GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,cellE,cellTime,cells);
1521  cellTime*=1e9;
1522  if(cellTime > 400)cellTime-=fConstantTimeShift;
1523 
1524  enerList[icell] = cellE;
1525  timeList[icell] = cellTime;
1526  }
1527 
1528  TMath::Sort(ncell,enerList,sortList);
1529 
1530  //
1531  // Location of max cell
1532  // First in sorted list.
1533  //
1534  Int_t absIdMax = cluster->GetCellAbsId(sortList[0]);
1535  // Float_t cellEMax = enerList[sortList[0]];
1536  Double_t cellTimeMax = timeList[sortList[0]];
1537 
1538  //
1539  // Cluster cell loop, select only secondary cells with enough contribution to cluster
1540  // Start from second highest energy cell
1541  //
1542  largeTime = 0;
1543  Float_t largeTimeE = 0;
1544  for(Int_t icell = 1; icell < ncell; icell++)
1545  {
1546  Int_t absId = cluster->GetCellAbsId(sortList[icell]);
1547  Float_t cellE = enerList[sortList[icell]];
1548  Double_t cellTime = timeList[sortList[icell]];
1549 
1550  //printf("\t icell %d, absId %d, ener %2.2f, time %2.2f, sort %d\n",icell,absId,cellE,cellTime,sortList[icell]);
1551 
1552  if ( absId == absIdMax ) continue;
1553 
1554  Float_t weight = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cellE,cluster->E());
1555 
1556  //printf("Cluster E %2.2f, cell E %2.2f, time %2.2f, weight %2.3f\n",cluster->E(),cellE,cellTime,weight);
1557  if(weight < 0.01) continue;
1558 
1559  if ( TMath::Abs(cellTime) > 50 )
1560  {
1561  largeTime++;
1562  largeTimeE += cellE;
1563  //printf(" cluster pT %2.2f, cell E %2.2f, large time %d; ",pt, cellE, largeTime);
1564  }
1565  if ( l0bin == -1 ) continue;
1566 
1567  fhTimeLam0BinPerSM [l0bin][sm]->Fill(pt, cellTime, GetEventWeight());
1568  fhTimeLam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellTime, weight*GetEventWeight());
1569 
1570  fhDTimeLam0BinPerSM [l0bin][sm]->Fill(pt, cellTimeMax-cellTime, GetEventWeight());
1571  fhDTimeLam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellTimeMax-cellTime, weight*GetEventWeight());
1572 
1573  fhCellClusterEFracLam0BinPerSM [l0bin][sm]->Fill(pt, cellE/cluster->E(), GetEventWeight());
1574 // fhCellClusterEFracLam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellE/cluster->E(), weight*GetEventWeight());
1575 
1576  fhCellClusterELam0BinPerSM [l0bin][sm]->Fill(pt, cellE, GetEventWeight());
1577  fhCellClusterELam0BinPerSMWeighted[l0bin][sm]->Fill(pt, cellE, weight*GetEventWeight());
1578 
1579  if(ptbin >= 0)
1580  {
1581  fhCellClusterIndexEAndTime[l0bin][ptbin] ->Fill(cellTime, icell , GetEventWeight());
1582  fhCellClusterEAndTime [l0bin][ptbin] ->Fill(cellTime, cellE , GetEventWeight());
1583  fhCellClusterEFracAndTime [l0bin][ptbin] ->Fill(cellTime, cellE/cluster->E(), GetEventWeight());
1584  }
1585 
1586  if ( TMath::Abs(cellTime) > 50 )
1587  {
1588  fhCellClusterEFracLam0BinPerSMLargeTime [l0bin][sm]->Fill(pt, cellE/cluster->E(), GetEventWeight());
1589  fhCellClusterELam0BinPerSMLargeTime [l0bin][sm]->Fill(pt, cellE , GetEventWeight());
1590  fhCellClusterIndexELam0BinPerSMLargeTime[l0bin][sm]->Fill(pt, icell , GetEventWeight());
1591  }
1592 // if(shared)
1593 // fhTimeLam0BinPerSMShared[l0bin][sm]->Fill(pt, cellTime, GetEventWeight());
1594 
1595  if(ptbin >= 0)
1596  {
1597  Int_t icol = -1, icolAbs = -1;
1598  Int_t irow = -1, irowAbs = -1;
1599  Int_t iRCU = -1;
1600  GetModuleNumberCellIndexesAbsCaloMap(absId,GetCalorimeter(), icol, irow, iRCU, icolAbs, irowAbs);
1601 
1602  fhColRowLam0BinPtBin [l0bin][ptbin]->Fill(icolAbs, irowAbs, GetEventWeight());
1603  fhColRowLam0BinPtBinWeighted[l0bin][ptbin]->Fill(icolAbs, irowAbs, weight*GetEventWeight());
1604 
1605  if ( TMath::Abs(cellTime) > 50 )
1606  fhColRowLam0BinPtBinLargeTime[l0bin][ptbin]->Fill(icolAbs, irowAbs, GetEventWeight());
1607 
1608 // if(shared) fhColRowLam0BinPtBinSMShared[l0bin][ptbin]->Fill(icolAbs, irowAbs, GetEventWeight());
1609  }
1610 
1611 
1612  } // cluster cell loop
1613  //if(largeTime > 0) printf("\n");
1614 
1615  if ( largeTime > 0 )
1616  {
1617  if ( l0bin != -1 )
1618  {
1619  if ( ptbin > 0 )
1620  fhEtaPhiLargeTimeInClusterCell[l0bin][ptbin]->Fill(eta, phi, GetEventWeight());
1621 
1622  fhCellClusterEFracLam0BinPerSMLargeTimeTotal[l0bin][sm]->Fill(pt, largeTimeE/cluster->E(), GetEventWeight());
1623  fhNCellsWithLargeTimeInCluster [l0bin][sm]->Fill(pt, largeTime , GetEventWeight());
1624  }
1625 
1626  fhLam0PerSMLargeTimeInClusterCell[sm]->Fill(pt, lambda0, GetEventWeight());
1627  fhLam1PerSMLargeTimeInClusterCell[sm]->Fill(pt, lambda1, GetEventWeight());
1628 
1629  if(largeTime < 6)
1630  {
1631  fhLam0PerNLargeTimeInClusterCell[largeTime-1]->Fill(pt, lambda0, GetEventWeight());
1632  fhLam1PerNLargeTimeInClusterCell[largeTime-1]->Fill(pt, lambda1, GetEventWeight());
1633  }
1634  }
1635 
1636  delete [] sortList ;
1637  delete [] enerList ;
1638  delete [] timeList ;
1639  //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);
1640  } // Check EMCal SM and regions
1641 
1642  Float_t l0 = 0., l1 = 0.;
1643  Float_t dispp= 0., dEta = 0., dPhi = 0.;
1644  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1646  {
1647  GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
1648  l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
1649  //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",
1650  // l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
1651  //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
1652  // disp, dPhi+dEta );
1653  fhDispEtaE -> Fill(energy, dEta , GetEventWeight());
1654  fhDispPhiE -> Fill(energy, dPhi , GetEventWeight());
1655  fhSumEtaE -> Fill(energy, sEta , GetEventWeight());
1656  fhSumPhiE -> Fill(energy, sPhi , GetEventWeight());
1657  fhSumEtaPhiE -> Fill(energy, sEtaPhi , GetEventWeight());
1658  fhDispEtaPhiDiffE -> Fill(energy, dPhi-dEta, GetEventWeight());
1659 
1660  if(dEta+dPhi>0)fhSphericityE -> Fill(energy,(dPhi-dEta)/(dEta+dPhi) , GetEventWeight());
1661  if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.), GetEventWeight());
1662  if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.), GetEventWeight());
1663 
1664  Int_t ebin = -1;
1665  if (energy < 2 ) ebin = 0;
1666  else if (energy < 4 ) ebin = 1;
1667  else if (energy < 6 ) ebin = 2;
1668  else if (energy < 10) ebin = 3;
1669  else if (energy < 15) ebin = 4;
1670  else if (energy < 20) ebin = 5;
1671  else ebin = 6;
1672 
1673  fhDispEtaDispPhi[ebin]->Fill(dEta , dPhi, GetEventWeight());
1674  fhLambda0DispEta[ebin]->Fill(lambda0, dEta, GetEventWeight());
1675  fhLambda0DispPhi[ebin]->Fill(lambda0, dPhi, GetEventWeight());
1676  }
1677 
1678  // If track-matching was off, check effect of track-matching residual cut
1679 
1680  if(!fRejectTrackMatch)
1681  {
1682  Float_t dEta = cluster->GetTrackDz();
1683  Float_t dPhi = cluster->GetTrackDx();
1684 // if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1685 // {
1686 // dPhi = 2000., dEta = 2000.;
1687 // GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dEta,dPhi);
1688 // }
1689 
1690  if(TMath::Abs(dEta) < 0.05 && TMath::Abs(dPhi) < 0.05)
1691  {
1692  fhLam0ETM ->Fill(energy, lambda0, GetEventWeight());
1693  fhLam0PtTM->Fill(pt , lambda0, GetEventWeight());
1694  fhLam1ETM ->Fill(energy, lambda1, GetEventWeight());
1696  fhDispETM ->Fill(energy, disp , GetEventWeight());
1697 
1698  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
1699  GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
1700  {
1701  fhLam0ETMTRD ->Fill(energy, lambda0, GetEventWeight());
1702  fhLam0PtTMTRD->Fill(pt , lambda0, GetEventWeight());
1703  fhLam1ETMTRD ->Fill(energy, lambda1, GetEventWeight());
1705  fhDispETMTRD ->Fill(energy, disp , GetEventWeight());
1706  }
1707  }
1708  } // If track-matching was off, check effect of matching residual cut
1709 
1711  {
1712  if(energy < 2)
1713  {
1714  fhNCellsLam0LowE ->Fill(ncells, lambda0, GetEventWeight());
1715  fhNCellsLam1LowE ->Fill(ncells, lambda1, GetEventWeight());
1716  fhNCellsDispLowE ->Fill(ncells, disp , GetEventWeight());
1717 
1718  fhLam1Lam0LowE ->Fill(lambda1, lambda0, GetEventWeight());
1719  fhLam0DispLowE ->Fill(lambda0, disp , GetEventWeight());
1720  fhDispLam1LowE ->Fill(disp , lambda1, GetEventWeight());
1721  fhEtaLam0LowE ->Fill(eta , lambda0, GetEventWeight());
1722  fhPhiLam0LowE ->Fill(phi , lambda0, GetEventWeight());
1723  }
1724  else
1725  {
1726  fhNCellsLam0HighE ->Fill(ncells, lambda0, GetEventWeight());
1727  fhNCellsLam1HighE ->Fill(ncells, lambda1, GetEventWeight());
1728  fhNCellsDispHighE ->Fill(ncells, disp , GetEventWeight());
1729 
1730  fhLam1Lam0HighE ->Fill(lambda1, lambda0, GetEventWeight());
1731  fhLam0DispHighE ->Fill(lambda0, disp , GetEventWeight());
1732  fhDispLam1HighE ->Fill(disp , lambda1, GetEventWeight());
1733  fhEtaLam0HighE ->Fill(eta , lambda0, GetEventWeight());
1734  fhPhiLam0HighE ->Fill(phi , lambda0, GetEventWeight());
1735  }
1736  }
1737 
1738  if(IsDataMC())
1739  {
1740  AliVCaloCells* cells = 0;
1741  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
1742  else cells = GetPHOSCells();
1743 
1744  // Fill histograms to check shape of embedded clusters
1745  Float_t fraction = 0;
1746  // printf("check embedding %i\n",GetReader()->IsEmbeddedClusterSelectionOn());
1747 
1748  // Compare the primary depositing more energy with the rest,
1749  // if no photon/electron as comon ancestor (conversions), count as other particle
1750  const UInt_t nlabels = cluster->GetNLabels();
1751  Int_t overpdg[nlabels];
1752  Int_t overlab[nlabels];
1753  Int_t noverlaps = 0;
1754  if(!GetReader()->IsEmbeddedClusterSelectionOn())
1755  noverlaps = GetMCAnalysisUtils()->GetNOverlaps(cluster->GetLabels(), nlabels,mcTag,-1,GetMC(),overpdg,overlab);
1756 
1757  //printf("N overlaps %d \n",noverlaps);
1758 
1759 
1761  {
1762  // Only working for EMCAL
1763  // printf("embedded\n");
1764 
1765  Float_t clusterE = 0; // recalculate in case corrections applied.
1766  Float_t cellE = 0;
1767  for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
1768  {
1769  cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
1770  clusterE+=cellE;
1771  fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
1772  }
1773 
1774  //Fraction of total energy due to the embedded signal
1775  fraction/=clusterE;
1776 
1777  AliDebug(1,Form("Energy fraction of embedded signal %2.3f, Energy %2.3f",fraction, clusterE));
1778 
1779  fhEmbeddedSignalFractionEnergy->Fill(clusterE, fraction, GetEventWeight());
1780  } // embedded fraction
1781 
1782  // Check the origin and fill histograms
1783 
1784  Int_t mcIndex = -1;
1785 
1786  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
1787  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
1788  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
1789  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1790  {
1791  mcIndex = kmcssPhoton ;
1792 
1793 // if(!GetReader()->IsEmbeddedClusterSelectionOn())
1794 // {
1795 // //Check particle overlaps in cluster
1796 //
1797 // if(noverlaps == 0)
1798 // {
1799 // fhMCPhotonELambda0NoOverlap ->Fill(energy, lambda0, GetEventWeight());
1800 // }
1801 // else if(noverlaps == 1)
1802 // {
1803 // fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0, GetEventWeight());
1804 // }
1805 // else if(noverlaps > 1)
1806 // {
1807 // fhMCPhotonELambda0NOverlap ->Fill(energy, lambda0, GetEventWeight());
1808 // }
1809 // else
1810 // {
1811 // AliWarning(Form("n overlaps = %d!!", noverlaps));
1812 // }
1813 // } // No embedding
1814 
1815  // Fill histograms to check shape of embedded clusters
1816  if(GetReader()->IsEmbeddedClusterSelectionOn())
1817  {
1818  if (fraction > 0.9)
1819  {
1820  fhEmbedPhotonELambda0FullSignal ->Fill(energy, lambda0, GetEventWeight());
1821  }
1822  else if(fraction > 0.5)
1823  {
1824  fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0, GetEventWeight());
1825  }
1826  else if(fraction > 0.1)
1827  {
1828  fhEmbedPhotonELambda0MostlyBkg ->Fill(energy, lambda0, GetEventWeight());
1829  }
1830  else
1831  {
1832  fhEmbedPhotonELambda0FullBkg ->Fill(energy, lambda0, GetEventWeight());
1833  }
1834  } // embedded
1835  } // photon no conversion
1836  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
1837  GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
1838  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
1839  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1840  {
1841  mcIndex = kmcssConversion ;
1842  }//conversion photon
1843 
1844  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
1845  {
1846  mcIndex = kmcssElectron ;
1847  }//electron
1848  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
1849  {
1850  mcIndex = kmcssPi0 ;
1851 
1852  // Fill histograms to check shape of embedded clusters
1853  if(GetReader()->IsEmbeddedClusterSelectionOn())
1854  {
1855  if (fraction > 0.9)
1856  {
1857  fhEmbedPi0ELambda0FullSignal ->Fill(energy, lambda0, GetEventWeight());
1858  }
1859  else if(fraction > 0.5)
1860  {
1861  fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0, GetEventWeight());
1862  }
1863  else if(fraction > 0.1)
1864  {
1865  fhEmbedPi0ELambda0MostlyBkg ->Fill(energy, lambda0, GetEventWeight());
1866  }
1867  else
1868  {
1869  fhEmbedPi0ELambda0FullBkg ->Fill(energy, lambda0, GetEventWeight());
1870  }
1871  } // embedded
1872 
1873  }//pi0
1874  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
1875  {
1876  mcIndex = kmcssEta ;
1877  }//eta
1878  else
1879  {
1880  mcIndex = kmcssOther ;
1881  }//other particles
1882 
1883  fhMCELambda0 [mcIndex]->Fill(energy, lambda0, GetEventWeight());
1884  fhMCPtLambda0 [mcIndex]->Fill(pt , lambda0, GetEventWeight());
1885  fhMCELambda1 [mcIndex]->Fill(energy, lambda1, GetEventWeight());
1886  fhMCNCellsE [mcIndex]->Fill(energy, ncells , GetEventWeight());
1887 
1889  {
1890  fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction, GetEventWeight());
1891  fhMCEDispersion [mcIndex]->Fill(energy, disp , GetEventWeight());
1892  }
1893 
1894  // Check particle overlaps in cluster
1895 
1896  if(!GetReader()->IsEmbeddedClusterSelectionOn())
1897  {
1898  if(noverlaps == 0)
1899  fhMCPtLambda0Overlaps[mcIndex][0]->Fill(pt, lambda0, GetEventWeight());
1900  else if(noverlaps == 1)
1901  fhMCPtLambda0Overlaps[mcIndex][1]->Fill(pt, lambda0, GetEventWeight());
1902  else if(noverlaps > 1)
1903  fhMCPtLambda0Overlaps[mcIndex][2]->Fill(pt, lambda0, GetEventWeight());
1904  else
1905  AliWarning(Form("n overlaps = %d!!", noverlaps));
1906  }
1907 
1909  {
1910  if (energy < 2.)
1911  {
1912  fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction, GetEventWeight());
1913  fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells, maxCellFraction, GetEventWeight());
1914  }
1915  else if(energy < 6.)
1916  {
1917  fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction, GetEventWeight());
1918  fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells, maxCellFraction, GetEventWeight());
1919  }
1920  else
1921  {
1922  fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction, GetEventWeight());
1923  fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells, maxCellFraction, GetEventWeight());
1924  }
1925 
1926  if(GetCalorimeter() == kEMCAL)
1927  {
1928  fhMCEDispEta [mcIndex]-> Fill(energy, dEta , GetEventWeight());
1929  fhMCEDispPhi [mcIndex]-> Fill(energy, dPhi , GetEventWeight());
1930  fhMCESumEtaPhi [mcIndex]-> Fill(energy, sEtaPhi , GetEventWeight());
1931  fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy, dPhi-dEta, GetEventWeight());
1932 
1933  if( dEta+dPhi > 0 ) fhMCESphericity[mcIndex]-> Fill(energy, (dPhi-dEta)/(dEta+dPhi), GetEventWeight());
1934 
1935  Int_t ebin = -1;
1936  if (energy < 2 ) ebin = 0;
1937  else if (energy < 4 ) ebin = 1;
1938  else if (energy < 6 ) ebin = 2;
1939  else if (energy < 10) ebin = 3;
1940  else if (energy < 15) ebin = 4;
1941  else if (energy < 20) ebin = 5;
1942  else ebin = 6;
1943 
1944  fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta , dPhi, GetEventWeight());
1945  fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0, dEta, GetEventWeight());
1946  fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0, dPhi, GetEventWeight());
1947  }
1948  }
1949  }// MC data
1950 }
1951 
1952 //__________________________________________________________________________
1959 //__________________________________________________________________________
1961  Int_t cut, Int_t tag)
1962 {
1963  Float_t dEta = cluster->GetTrackDz();
1964  Float_t dPhi = cluster->GetTrackDx();
1965 
1966 // if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1967 // {
1968 // dPhi = 2000., dEta = 2000.;
1969 // GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedPhiesiduals(cluster->GetID(),dEta,dPhi);
1970 // }
1971 
1972  AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1973 
1974  Bool_t positive = kFALSE;
1975  if(track) positive = (track->Charge()>0);
1976 
1977  if(fhTrackMatchedDEtaPos[cut] && TMath::Abs(dPhi) < 999)
1978  {
1979 // fhTrackMatchedDEta[cut]->Fill(cluster->E(), dEta, GetEventWeight());
1980 // fhTrackMatchedDPhi[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
1981 // if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dEta, dPhi, GetEventWeight());
1982 
1983  if(track)
1984  {
1985 // fhTrackMatchedDEtaTrackPt[cut]->Fill(track->Pt(), dEta, GetEventWeight());
1986 // fhTrackMatchedDPhiTrackPt[cut]->Fill(track->Pt(), dPhi, GetEventWeight());
1987 // if(track->Pt() > 0.5) fhTrackMatchedDEtaDPhiTrackPt[cut]->Fill(dEta, dPhi, GetEventWeight());
1988 
1989  if(positive)
1990  {
1991  fhTrackMatchedDEtaPos[cut]->Fill(cluster->E(), dEta, GetEventWeight());
1992  fhTrackMatchedDPhiPos[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
1993  if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiPos[cut]->Fill(dEta,dPhi, GetEventWeight());
1994 
1996  {
1997  fhTrackMatchedDEtaPosTrackPt[cut]->Fill(track->Pt(), dEta, GetEventWeight());
1998  fhTrackMatchedDPhiPosTrackPt[cut]->Fill(track->Pt(), dPhi, GetEventWeight());
1999  if(track->Pt() > 0.5) fhTrackMatchedDEtaDPhiPosTrackPt[cut]->Fill(dEta,dPhi, GetEventWeight());
2000  }
2001  }
2002  else
2003  {
2004  fhTrackMatchedDEtaNeg[cut]->Fill(cluster->E(), dEta, GetEventWeight());
2005  fhTrackMatchedDPhiNeg[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
2006  if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiNeg[cut]->Fill(dEta, dPhi, GetEventWeight());
2007 
2009  {
2010  fhTrackMatchedDEtaNegTrackPt[cut]->Fill(track->Pt(), dEta, GetEventWeight());
2011  fhTrackMatchedDPhiNegTrackPt[cut]->Fill(track->Pt(), dPhi, GetEventWeight());
2012  if(track->Pt() > 0.5) fhTrackMatchedDEtaDPhiNegTrackPt[cut]->Fill(dEta, dPhi, GetEventWeight());
2013  }
2014  }
2015  }
2016 
2017  Int_t nSMod = GetModuleNumber(cluster);
2018 
2019  if( GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
2020  nSMod >= GetFirstSMCoveredByTRD() )
2021  {
2022  fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(), dEta, GetEventWeight());
2023  fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
2024  }
2025 
2026  // Check dEdx and E/p of matched clusters
2027 
2028  if(TMath::Abs(dEta) < 0.05 && TMath::Abs(dPhi) < 0.05)
2029  {
2030  if(track)
2031  {
2032  Float_t dEdx = track->GetTPCsignal();
2033  Float_t eOverp = cluster->E()/track->P();
2034 
2035  fhdEdx [cut]->Fill(cluster->E(), dEdx , GetEventWeight());
2036  fhEOverP[cut]->Fill(cluster->E(), eOverp, GetEventWeight());
2037 
2039  {
2040  fhdEdxTrackPt [cut]->Fill(track->Pt(), dEdx , GetEventWeight());
2041  fhEOverPTrackPt[cut]->Fill(track->Pt(), eOverp, GetEventWeight());
2042  }
2043 
2044  if(GetCalorimeter()==kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
2045  nSMod >= GetFirstSMCoveredByTRD() )
2046  fhEOverPTRD[cut]->Fill(cluster->E(), eOverp, GetEventWeight());
2047  }
2048  else
2049  AliWarning(Form("Residual OK but (dPhi, dEta)= (%2.4f,%2.4f) no track associated WHAT?", dPhi,dEta));
2050 
2051  if(IsDataMC())
2052  {
2053  if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2054  {
2055  if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
2056  GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
2057  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5, GetEventWeight());
2059  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5, GetEventWeight());
2061  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5, GetEventWeight());
2062  else
2063  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5, GetEventWeight());
2064 
2065  // Check if several particles contributed to cluster and discard overlapped mesons
2068  {
2069  if(cluster->GetNLabels()==1)
2070  {
2071  fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(), dEta, GetEventWeight());
2072  fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
2073  }
2074  else
2075  {
2076  fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(), dEta, GetEventWeight());
2077  fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
2078  }
2079 
2080  }// Check overlaps
2081 
2082  }
2083  else
2084  {
2087  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5, GetEventWeight());
2089  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5, GetEventWeight());
2091  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5, GetEventWeight());
2092  else
2093  fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5, GetEventWeight());
2094 
2095  // Check if several particles contributed to cluster
2098  {
2099  fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(), dEta, GetEventWeight());
2100  fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(), dPhi, GetEventWeight());
2101 
2102  }// Check overlaps
2103 
2104  }
2105 
2106  } // MC
2107 
2108  } // residuals window
2109 
2110  } // Small residual
2111 }
2112 
2113 //___________________________________________
2115 //___________________________________________
2117 {
2118  TString parList ; //this will be list of parameters used for this analysis.
2119  const Int_t buffersize = 255;
2120  char onePar[buffersize] ;
2121 
2122  snprintf(onePar,buffersize,"--- AliAnaPhoton ---:") ;
2123  parList+=onePar ;
2124  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
2125  parList+=onePar ;
2126  snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster);",fMinDist) ;
2127  parList+=onePar ;
2128  snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation);",fMinDist2) ;
2129  parList+=onePar ;
2130  snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study);",fMinDist3) ;
2131  parList+=onePar ;
2132  snprintf(onePar,buffersize,"fRejectTrackMatch: %d",fRejectTrackMatch) ;
2133  parList+=onePar ;
2134 
2135  // Get parameters set in base class.
2136  parList += GetBaseParametersList() ;
2137 
2138  // Get parameters set in PID class.
2139  parList += GetCaloPID()->GetPIDParametersList() ;
2140 
2141  // Get parameters set in FiducialCut class (not available yet)
2142  //parlist += GetFidCut()->GetFidCutParametersList()
2143 
2144  return new TObjString(parList) ;
2145 }
2146 
2147 //_____________________________________________
2150 //_____________________________________________
2152 {
2153  TList * outputContainer = new TList() ;
2154  outputContainer->SetName("PhotonHistos") ;
2155 
2162 
2169 
2170  Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
2173  Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
2174  Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
2175  Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
2176 
2177  Int_t nratbins = GetHistogramRanges()->GetHistoRatioBins();
2180  Int_t ndifbins = GetHistogramRanges()->GetHistoEDiffBins();
2183 
2187 
2191 
2192  // Init the number of modules, set in the class AliCalorimeterUtils
2193  //
2194  InitCaloParameters(); // See AliCaloTrackCorrBaseClass
2195 
2196  Int_t totalSM = fLastModule-fFirstModule+1;
2197 
2198  //printf("N SM %d, first SM %d, last SM %d, total %d\n",fNModules,fFirstModule,fLastModule, totalSM);
2199 
2200  // Cell column-row histograms, see base class for data members setting
2201  //fNMaxColsFull+2,-1.5,fNMaxColsFull+0.5, fNMaxRowsFull+2,-1.5,fNMaxRowsFull+0.5
2202  Int_t ncolcell = fNMaxColsFull+2;
2203  Float_t colcellmin = -1.5;
2204  Float_t colcellmax = fNMaxColsFull+0.5;
2205 
2207  Float_t rowcellmin = fNMaxRowsFullMin-1.5;
2208  Float_t rowcellmax = fNMaxRowsFullMax+0.5;
2209 
2210  Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins for SS studies
2211 
2212  TString cut[] = {"Open","Reader","E","Time","NCells","NLM","Fidutial","Matching","Bad","PID"};
2213  for (Int_t i = 0; i < 10 ; i++)
2214  {
2215  fhClusterCutsE[i] = new TH1F(Form("hE_Cut_%d_%s", i, cut[i].Data()),
2216  Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
2217  nptbins,ptmin,ptmax);
2218  fhClusterCutsE[i]->SetYTitle("d#it{N}/d#it{E} ");
2219  fhClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
2220  outputContainer->Add(fhClusterCutsE[i]) ;
2221 
2222  fhClusterCutsPt[i] = new TH1F(Form("hPt_Cut_%d_%s", i, cut[i].Data()),
2223  Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
2224  nptbins,ptmin,ptmax);
2225  fhClusterCutsPt[i]->SetYTitle("d#it{N}/d#it{E} ");
2226  fhClusterCutsPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2227  outputContainer->Add(fhClusterCutsPt[i]) ;
2228  }
2229 
2230  fhEClusterSM = new TH2F("hEClusterSM","Raw clusters E and super-module number",
2231  nptbins,ptmin,ptmax,
2232  totalSM,fFirstModule-0.5,fLastModule+0.5);
2233  fhEClusterSM->SetYTitle("SuperModule ");
2234  fhEClusterSM->SetXTitle("#it{E} (GeV)");
2235  outputContainer->Add(fhEClusterSM) ;
2236 
2237  fhPtClusterSM = new TH2F("hPtClusterSM","Raw clusters #it{p}_{T} and super-module number",
2238  nptbins,ptmin,ptmax,
2239  totalSM,fFirstModule-0.5,fLastModule+0.5);
2240  fhPtClusterSM->SetYTitle("SuperModule ");
2241  fhPtClusterSM->SetXTitle("#it{E} (GeV)");
2242  outputContainer->Add(fhPtClusterSM) ;
2243 
2244  fhEPhotonSM = new TH2F("hEPhotonSM","Selected clusters E and super-module number",
2245  nptbins,ptmin,ptmax,
2246  totalSM,fFirstModule-0.5,fLastModule+0.5);
2247  fhEPhotonSM->SetYTitle("SuperModule ");
2248  fhEPhotonSM->SetXTitle("#it{E} (GeV)");
2249  outputContainer->Add(fhEPhotonSM) ;
2250 
2251  fhPtPhotonSM = new TH2F("hPtPhotonSM","Selected clusters #it{p}_{T} and super-module number",
2252  nptbins,ptmin,ptmax,
2253  totalSM,fFirstModule-0.5,fLastModule+0.5);
2254  fhPtPhotonSM->SetYTitle("SuperModule ");
2255  fhPtPhotonSM->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2256  outputContainer->Add(fhPtPhotonSM) ;
2257 
2259  {
2260 // printf("Bins Mult: n %d, min %d, max %d; Sum: n %d, min %f, max %f\n",
2261 // nmultbin,multmin,multmax,nsumbin,summin, summax);
2262  for(Int_t icut = 0; icut < GetReader()->GetTrackMultiplicityNPtCut(); icut++)
2263  {
2264  fhPtPhotonNTracks[icut] = new TH2F
2265  (Form("hPtPhotonNTracks_PtCut%d",icut),
2266  Form("Number of tracks per event with |#eta|<%2.2f and #it{p}_{T} > %2.2f GeV/#it{c}",
2267  GetReader()->GetTrackMultiplicityEtaCut(),GetReader()->GetTrackMultiplicityPtCut(icut)),
2268  nptbins,ptmin,ptmax, nmultbin,multmin,multmax);
2269  fhPtPhotonNTracks[icut]->SetYTitle("# of tracks");
2270  fhPtPhotonNTracks[icut]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2271  outputContainer->Add(fhPtPhotonNTracks[icut]);
2272 
2273  fhPtPhotonSumPtTracks[icut] = new TH2F
2274  (Form("hPtPhotonSumPtTracks_PtCut%d",icut),
2275  Form("#Sigma #it{p}_{T} of tracks per event with |#eta|<%2.2f and #it{p}_{T} > %2.2f GeV/#it{c}",
2276  GetReader()->GetTrackMultiplicityEtaCut(),GetReader()->GetTrackMultiplicityPtCut(icut)),
2277  nptbins,ptmin,ptmax, nsumbin,summin,summax);
2278  fhPtPhotonSumPtTracks[icut]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2279  fhPtPhotonSumPtTracks[icut]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2280  outputContainer->Add(fhPtPhotonSumPtTracks[icut]);
2281  }
2282  }
2283 
2284  fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
2285  fhNCellsE->SetXTitle("#it{E} (GeV)");
2286  fhNCellsE->SetYTitle("# of cells in cluster");
2287  outputContainer->Add(fhNCellsE);
2288 
2289  fhCellsE = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
2290  fhCellsE->SetXTitle("#it{E}_{cluster} (GeV)");
2291  fhCellsE->SetYTitle("#it{E}_{cell} (GeV)");
2292  outputContainer->Add(fhCellsE);
2293 
2294  fhTimePt = new TH2F ("hTimePt","time of cluster vs pT of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2295  fhTimePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2296  fhTimePt->SetYTitle("#it{time} (ns)");
2297  outputContainer->Add(fhTimePt);
2298 
2300  {
2301  fhMaxCellDiffClusterE = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
2302  nptbins,ptmin,ptmax, 500,0,1.);
2303  fhMaxCellDiffClusterE->SetXTitle("#it{E}_{cluster} (GeV) ");
2304  fhMaxCellDiffClusterE->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
2305  outputContainer->Add(fhMaxCellDiffClusterE);
2306  }
2307 
2308  fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
2309  fhEPhoton->SetYTitle("#it{counts}");
2310  fhEPhoton->SetXTitle("#it{E}_{#gamma}(GeV)");
2311  outputContainer->Add(fhEPhoton) ;
2312 
2313  fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs #it{p}_{T}",nptbins,ptmin,ptmax);
2314  fhPtPhoton->SetYTitle("#it{counts}");
2315  fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/#it{c})");
2316  outputContainer->Add(fhPtPhoton) ;
2317 
2319  {
2320  fhPtCentralityPhoton = new TH2F("hPtCentralityPhoton","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
2321  fhPtCentralityPhoton->SetYTitle("Centrality");
2322  fhPtCentralityPhoton->SetXTitle("#it{p}_{T}(GeV/#it{c})");
2323  outputContainer->Add(fhPtCentralityPhoton) ;
2324 
2325  fhPtEventPlanePhoton = new TH2F("hPtEventPlanePhoton","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
2326  fhPtEventPlanePhoton->SetYTitle("Event plane angle (rad)");
2327  fhPtEventPlanePhoton->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2328  outputContainer->Add(fhPtEventPlanePhoton) ;
2329  }
2330 
2331  fhEtaPhi = new TH2F
2332  ("hEtaPhi","cluster,#it{E} > 0.5 GeV, #eta vs #varphi",netabins,etamin,etamax,nphibins,phimin,phimax);
2333  fhEtaPhi->SetYTitle("#varphi (rad)");
2334  fhEtaPhi->SetXTitle("#eta");
2335  outputContainer->Add(fhEtaPhi) ;
2336 
2337  fhPhiPhoton = new TH2F
2338  ("hPhiPhoton","#varphi_{#gamma} vs #it{p}_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2339  fhPhiPhoton->SetYTitle("#varphi (rad)");
2340  fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
2341  outputContainer->Add(fhPhiPhoton) ;
2342 
2343  fhEtaPhoton = new TH2F
2344  ("hEtaPhoton","#eta_{#gamma} vs #it{p}_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2345  fhEtaPhoton->SetYTitle("#eta");
2346  fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
2347  outputContainer->Add(fhEtaPhoton) ;
2348 
2349  fhEtaPhiPhoton = new TH2F
2350  ("hEtaPhiPhoton","#eta vs #varphi",netabins,etamin,etamax,nphibins,phimin,phimax);
2351  fhEtaPhiPhoton->SetYTitle("#varphi (rad)");
2352  fhEtaPhiPhoton->SetXTitle("#eta");
2353  outputContainer->Add(fhEtaPhiPhoton) ;
2354  if(GetMinPt() < 0.5)
2355  {
2356  fhEtaPhi05Photon = new TH2F
2357  ("hEtaPhi05Photon","#eta vs #varphi, E < 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
2358  fhEtaPhi05Photon->SetYTitle("#varphi (rad)");
2359  fhEtaPhi05Photon->SetXTitle("#eta");
2360  outputContainer->Add(fhEtaPhi05Photon) ;
2361  }
2362 
2363  fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
2364  nptbins,ptmin,ptmax,10,0,10);
2365  fhNLocMax ->SetYTitle("N maxima");
2366  fhNLocMax ->SetXTitle("#it{E} (GeV)");
2367  outputContainer->Add(fhNLocMax) ;
2368 
2369  //Shower shape
2370  if(fFillSSHistograms)
2371  {
2372  fhLam0E = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2373  fhLam0E->SetYTitle("#lambda_{0}^{2}");
2374  fhLam0E->SetXTitle("#it{E} (GeV)");
2375  outputContainer->Add(fhLam0E);
2376 
2377  fhLam0Pt = new TH2F ("hLam0Pt","#lambda_{0}^{2} vs #it{p}_{T}", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2378  fhLam0Pt->SetYTitle("#lambda_{0}^{2}");
2379  fhLam0Pt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2380  outputContainer->Add(fhLam0Pt);
2381 
2382  fhLam1E = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2383  fhLam1E->SetYTitle("#lambda_{1}^{2}");
2384  fhLam1E->SetXTitle("#it{E} (GeV)");
2385  outputContainer->Add(fhLam1E);
2386 
2387  fhLam1Pt = new TH2F ("hLam1Pt","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2388  fhLam1Pt->SetYTitle("#lambda_{1}^{2}");
2389  fhLam1Pt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2390  outputContainer->Add(fhLam1Pt);
2391 
2393  {
2394  fhDispE = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2395  fhDispE->SetYTitle("D^{2}");
2396  fhDispE->SetXTitle("#it{E} (GeV) ");
2397  outputContainer->Add(fhDispE);
2398 
2399  fhDispPt = new TH2F ("hDispPt"," dispersion^{2} vs #it{p}_{T}", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2400  fhDispPt->SetYTitle("D^{2}");
2401  fhDispPt->SetXTitle("#it{p}_{T} (GeV/#it{c}) ");
2402  outputContainer->Add(fhDispPt);
2403  }
2404 
2406  {
2407  fhLam0PtNLM1 = new TH2F ("hLam0PtNLM1","#lambda_{0}^{2} vs #it{p}_{T}, #it{n}_{LM}=1", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2408  fhLam0PtNLM1->SetYTitle("#lambda_{0}^{2}");
2409  fhLam0PtNLM1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2410  outputContainer->Add(fhLam0PtNLM1);
2411 
2412  fhLam0PtNLM2 = new TH2F ("hLam0PtNLM2","#lambda_{0}^{2} vs #it{p}_{T}, #it{n}_{LM}=2", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2413  fhLam0PtNLM2->SetYTitle("#lambda_{0}^{2}");
2414  fhLam0PtNLM2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2415  outputContainer->Add(fhLam0PtNLM2);
2416 
2417  fhLam1PtNLM1 = new TH2F ("hLam1PtNLM1","#lambda_{1}^{2} vs #it{p}_{T}, #it{n}_{LM}=1", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2418  fhLam1PtNLM1->SetYTitle("#lambda_{1}^{2}");
2419  fhLam1PtNLM1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2420  outputContainer->Add(fhLam1PtNLM1);
2421 
2422  fhLam1PtNLM2 = new TH2F ("hLam1PtNLM2","#lambda_{1}^{2} vs #it{p}_{T}, #it{n}_{LM}=2", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2423  fhLam1PtNLM2->SetYTitle("#lambda_{1}^{2}");
2424  fhLam1PtNLM2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2425  outputContainer->Add(fhLam1PtNLM2);
2426  }
2427 
2428  if(!fRejectTrackMatch)
2429  {
2430  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);
2431  fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
2432  fhLam0ETM->SetXTitle("#it{E} (GeV)");
2433  outputContainer->Add(fhLam0ETM);
2434 
2435  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);
2436  fhLam0PtTM->SetYTitle("#lambda_{0}^{2}");
2437  fhLam0PtTM->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2438  outputContainer->Add(fhLam0PtTM);
2439 
2440  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);
2441  fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
2442  fhLam1ETM->SetXTitle("#it{E} (GeV)");
2443  outputContainer->Add(fhLam1ETM);
2444 
2446  {
2447  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);
2448  fhDispETM->SetYTitle("D^{2}");
2449  fhDispETM->SetXTitle("#it{E} (GeV) ");
2450  outputContainer->Add(fhDispETM);
2451  }
2452  }
2453 
2454  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
2455  {
2456  fhLam0ETRD = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2457  fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
2458  fhLam0ETRD->SetXTitle("#it{E} (GeV)");
2459  outputContainer->Add(fhLam0ETRD);
2460 
2461  fhLam0PtTRD = new TH2F ("hLam0PtTRD","#lambda_{0}^{2} vs #it{p}_{T}, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2462  fhLam0PtTRD->SetYTitle("#lambda_{0}^{2}");
2463  fhLam0PtTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2464  outputContainer->Add(fhLam0PtTRD);
2465 
2466  fhLam1ETRD = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2467  fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
2468  fhLam1ETRD->SetXTitle("#it{E} (GeV)");
2469  outputContainer->Add(fhLam1ETRD);
2470 
2472  {
2473  fhDispETRD = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2474  fhDispETRD->SetYTitle("Dispersion^{2}");
2475  fhDispETRD->SetXTitle("#it{E} (GeV) ");
2476  outputContainer->Add(fhDispETRD);
2477  }
2478 
2480  {
2481  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);
2482  fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
2483  fhLam0ETMTRD->SetXTitle("#it{E} (GeV)");
2484  outputContainer->Add(fhLam0ETMTRD);
2485 
2486  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);
2487  fhLam0PtTMTRD->SetYTitle("#lambda_{0}^{2}");
2488  fhLam0PtTMTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2489  outputContainer->Add(fhLam0PtTMTRD);
2490 
2491  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);
2492  fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
2493  fhLam1ETMTRD->SetXTitle("#it{E} (GeV)");
2494  outputContainer->Add(fhLam1ETMTRD);
2495 
2497  {
2498  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);
2499  fhDispETMTRD->SetYTitle("Dispersion^{2}");
2500  fhDispETMTRD->SetXTitle("#it{E} (GeV) ");
2501  outputContainer->Add(fhDispETMTRD);
2502  }
2503  }
2504  }
2505 
2507  {
2508  fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2509  fhNCellsLam0LowE->SetXTitle("N_{cells}");
2510  fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
2511  outputContainer->Add(fhNCellsLam0LowE);
2512 
2513  fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2514  fhNCellsLam0HighE->SetXTitle("N_{cells}");
2515  fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
2516  outputContainer->Add(fhNCellsLam0HighE);
2517 
2518  fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2519  fhNCellsLam1LowE->SetXTitle("N_{cells}");
2520  fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
2521  outputContainer->Add(fhNCellsLam1LowE);
2522 
2523  fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2524  fhNCellsLam1HighE->SetXTitle("N_{cells}");
2525  fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
2526  outputContainer->Add(fhNCellsLam1HighE);
2527 
2528  fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2529  fhNCellsDispLowE->SetXTitle("N_{cells}");
2530  fhNCellsDispLowE->SetYTitle("D^{2}");
2531  outputContainer->Add(fhNCellsDispLowE);
2532 
2533  fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2534  fhNCellsDispHighE->SetXTitle("N_{cells}");
2535  fhNCellsDispHighE->SetYTitle("D^{2}");
2536  outputContainer->Add(fhNCellsDispHighE);
2537 
2538  fhEtaLam0LowE = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
2539  fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
2540  fhEtaLam0LowE->SetXTitle("#eta");
2541  outputContainer->Add(fhEtaLam0LowE);
2542 
2543  fhPhiLam0LowE = new TH2F ("hPhiLam0LowE","#varphi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
2544  fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
2545  fhPhiLam0LowE->SetXTitle("#varphi (rad)");
2546  outputContainer->Add(fhPhiLam0LowE);
2547 
2548  fhEtaLam0HighE = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, #it{E} > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
2549  fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
2550  fhEtaLam0HighE->SetXTitle("#eta");
2551  outputContainer->Add(fhEtaLam0HighE);
2552 
2553  fhPhiLam0HighE = new TH2F ("hPhiLam0HighE","#varphi vs #lambda_{0}^{2}, #it{E} > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
2554  fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
2555  fhPhiLam0HighE->SetXTitle("#varphi (rad)");
2556  outputContainer->Add(fhPhiLam0HighE);
2557 
2558  fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2559  fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
2560  fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
2561  outputContainer->Add(fhLam1Lam0LowE);
2562 
2563  fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2564  fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
2565  fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
2566  outputContainer->Add(fhLam1Lam0HighE);
2567 
2568  fhLam0DispLowE = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2569  fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
2570  fhLam0DispLowE->SetYTitle("D^{2}");
2571  outputContainer->Add(fhLam0DispLowE);
2572 
2573  fhLam0DispHighE = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2574  fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
2575  fhLam0DispHighE->SetYTitle("D^{2}");
2576  outputContainer->Add(fhLam0DispHighE);
2577 
2578  fhDispLam1LowE = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2579  fhDispLam1LowE->SetXTitle("D^{2}");
2580  fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
2581  outputContainer->Add(fhDispLam1LowE);
2582 
2583  fhDispLam1HighE = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2584  fhDispLam1HighE->SetXTitle("D^{2}");
2585  fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
2586  outputContainer->Add(fhDispLam1HighE);
2587 
2588  if(GetCalorimeter() == kEMCAL)
2589  {
2590  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);
2591  fhDispEtaE->SetXTitle("#it{E} (GeV)");
2592  fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
2593  outputContainer->Add(fhDispEtaE);
2594 
2595  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);
2596  fhDispPhiE->SetXTitle("#it{E} (GeV)");
2597  fhDispPhiE->SetYTitle("#sigma^{2}_{#varphi #varphi}");
2598  outputContainer->Add(fhDispPhiE);
2599 
2600  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);
2601  fhSumEtaE->SetXTitle("#it{E} (GeV)");
2602  fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
2603  outputContainer->Add(fhSumEtaE);
2604 
2605  fhSumPhiE = new TH2F ("hSumPhiE","#delta^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i})^{2}/ #Sigma w_{i} - <#varphi>^{2} vs E",
2606  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2607  fhSumPhiE->SetXTitle("#it{E} (GeV)");
2608  fhSumPhiE->SetYTitle("#delta^{2}_{#varphi #varphi}");
2609  outputContainer->Add(fhSumPhiE);
2610 
2611  fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #varphi} = #Sigma w_{i}(#varphi_{i} #eta_{i} ) / #Sigma w_{i} - <#varphi><#eta> vs E",
2612  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
2613  fhSumEtaPhiE->SetXTitle("#it{E} (GeV)");
2614  fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #varphi}");
2615  outputContainer->Add(fhSumEtaPhiE);
2616 
2617  fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta} vs E",
2618  nptbins,ptmin,ptmax,200, -10,10);
2619  fhDispEtaPhiDiffE->SetXTitle("#it{E} (GeV)");
2620  fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#varphi #varphi}-#sigma^{2}_{#eta #eta}");
2621  outputContainer->Add(fhDispEtaPhiDiffE);
2622 
2623  fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi}) vs E",
2624  nptbins,ptmin,ptmax, 200, -1,1);
2625  fhSphericityE->SetXTitle("#it{E} (GeV)");
2626  fhSphericityE->SetYTitle("s = (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi})");
2627  outputContainer->Add(fhSphericityE);
2628 
2629  fhDispSumEtaDiffE = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
2630  fhDispSumEtaDiffE->SetXTitle("#it{E} (GeV)");
2631  fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average");
2632  outputContainer->Add(fhDispSumEtaDiffE);
2633 
2634  fhDispSumPhiDiffE = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#varphi #varphi} - #delta^{2}_{#varphi #varphi} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
2635  fhDispSumPhiDiffE->SetXTitle("#it{E} (GeV)");
2636  fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#varphi #varphi} - #delta^{2}_{#varphi #varphi} / average");
2637  outputContainer->Add(fhDispSumPhiDiffE);
2638 
2639  for(Int_t i = 0; i < 7; i++)
2640  {
2641  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]),
2642  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2643  fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2644  fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
2645  outputContainer->Add(fhDispEtaDispPhi[i]);
2646 
2647  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]),
2648  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2649  fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
2650  fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
2651  outputContainer->Add(fhLambda0DispEta[i]);
2652 
2653  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]),
2654  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2655  fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
2656  fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
2657  outputContainer->Add(fhLambda0DispPhi[i]);
2658  }
2659  }
2660  }
2661  } // Shower shape
2662 
2663  // Track Matching
2664 
2665  if(fFillTMHisto)
2666  {
2667  TString cutTM [] = {"NoCut",""};
2668 
2669  for(Int_t i = 0; i < 2; i++)
2670  {
2671 // fhTrackMatchedDEta[i] = new TH2F
2672 // (Form("hTrackMatchedDEta%s",cutTM[i].Data()),
2673 // Form("d#eta of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2674 // nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2675 // fhTrackMatchedDEta[i]->SetYTitle("#Delta #eta");
2676 // fhTrackMatchedDEta[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2677 //
2678 // fhTrackMatchedDPhi[i] = new TH2F
2679 // (Form("hTrackMatchedDPhi%s",cutTM[i].Data()),
2680 // Form("#Delta #varphi of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2681 // nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2682 // fhTrackMatchedDPhi[i]->SetYTitle("#Delta #varphi (rad)");
2683 // fhTrackMatchedDPhi[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2684 //
2685 // fhTrackMatchedDEtaDPhi[i] = new TH2F
2686 // (Form("hTrackMatchedDEtaDPhi%s",cutTM[i].Data()),
2687 // Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{E}_{cluster}, %s",cutTM[i].Data()),
2688 // nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2689 // fhTrackMatchedDEtaDPhi[i]->SetYTitle("#Delta #varphi (rad)");
2690 // fhTrackMatchedDEtaDPhi[i]->SetXTitle("#Delta #eta");
2691 
2692  fhTrackMatchedDEtaPos[i] = new TH2F
2693  (Form("hTrackMatchedDEtaPos%s",cutTM[i].Data()),
2694  Form("#Delta #eta of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2695  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2696  fhTrackMatchedDEtaPos[i]->SetYTitle("#Delta #eta");
2697  fhTrackMatchedDEtaPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2698 
2699  fhTrackMatchedDPhiPos[i] = new TH2F
2700  (Form("hTrackMatchedDPhiPos%s",cutTM[i].Data()),
2701  Form("#Delta #varphi of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2702  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2703  fhTrackMatchedDPhiPos[i]->SetYTitle("#Delta #varphi (rad)");
2704  fhTrackMatchedDPhiPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2705 
2707  (Form("hTrackMatchedDEtaDPhiPos%s",cutTM[i].Data()),
2708  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{E}_{cluster}, %s",cutTM[i].Data()),
2709  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2710  fhTrackMatchedDEtaDPhiPos[i]->SetYTitle("#Delta #varphi (rad)");
2711  fhTrackMatchedDEtaDPhiPos[i]->SetXTitle("#Delta #eta");
2712 
2713  fhTrackMatchedDEtaNeg[i] = new TH2F
2714  (Form("hTrackMatchedDEtaNeg%s",cutTM[i].Data()),
2715  Form("#Delta #eta of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2716  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2717  fhTrackMatchedDEtaNeg[i]->SetYTitle("#Delta #eta");
2718  fhTrackMatchedDEtaNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2719 
2720  fhTrackMatchedDPhiNeg[i] = new TH2F
2721  (Form("hTrackMatchedDPhiNeg%s",cutTM[i].Data()),
2722  Form("#Delta #varphi of cluster-track vs #it{E}_{cluster}, %s",cutTM[i].Data()),
2723  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2724  fhTrackMatchedDPhiNeg[i]->SetYTitle("#Delta #varphi (rad)");
2725  fhTrackMatchedDPhiNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2726 
2728  (Form("hTrackMatchedDEtaDPhiNeg%s",cutTM[i].Data()),
2729  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{E}_{cluster} > 0.5 GeV, %s",cutTM[i].Data()),
2730  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2731  fhTrackMatchedDEtaDPhiNeg[i]->SetYTitle("#Delta #varphi (rad)");
2732  fhTrackMatchedDEtaDPhiNeg[i]->SetXTitle("#Delta #eta");
2733 
2734  fhdEdx[i] = new TH2F (Form("hdEdx%s",cutTM[i].Data()),Form("matched track <dE/dx> vs cluster #it{E}, %s",cutTM[i].Data()),
2735  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2736  fhdEdx[i]->SetXTitle("#it{E}^{cluster} (GeV)");
2737  fhdEdx[i]->SetYTitle("<dE/dx>");
2738 
2739  fhEOverP[i] = new TH2F (Form("hEOverP%s",cutTM[i].Data()),Form("matched track E/p vs cluster #it{E}, %s",cutTM[i].Data()),
2740  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2741  fhEOverP[i]->SetXTitle("#it{E}^{cluster} (GeV)");
2742  fhEOverP[i]->SetYTitle("E/p");
2743 
2744 // outputContainer->Add(fhTrackMatchedDEta[i]) ;
2745 // outputContainer->Add(fhTrackMatchedDPhi[i]) ;
2746 // outputContainer->Add(fhTrackMatchedDEtaDPhi[i]) ;
2747  outputContainer->Add(fhTrackMatchedDEtaPos[i]) ;
2748  outputContainer->Add(fhTrackMatchedDPhiPos[i]) ;
2749  outputContainer->Add(fhTrackMatchedDEtaDPhiPos[i]) ;
2750  outputContainer->Add(fhTrackMatchedDEtaNeg[i]) ;
2751  outputContainer->Add(fhTrackMatchedDPhiNeg[i]) ;
2752  outputContainer->Add(fhTrackMatchedDEtaDPhiNeg[i]) ;
2753  outputContainer->Add(fhdEdx[i]);
2754  outputContainer->Add(fhEOverP[i]);
2755 
2757  {
2758 // fhTrackMatchedDEtaTrackPt[i] = new TH2F
2759 // (Form("hTrackMatchedDEtaTrackPt%s",cutTM[i].Data()),
2760 // Form("#Delta #eta of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2761 // nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2762 // fhTrackMatchedDEtaTrackPt[i]->SetYTitle("#Delta #eta");
2763 // fhTrackMatchedDEtaTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2764 //
2765 // fhTrackMatchedDPhiTrackPt[i] = new TH2F
2766 // (Form("hTrackMatchedDPhiTrackPt%s",cutTM[i].Data()),
2767 // Form("#Delta #varphi of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2768 // nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2769 // fhTrackMatchedDPhiTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2770 // fhTrackMatchedDPhiTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2771 //
2772 // fhTrackMatchedDEtaDPhiTrackPt[i] = new TH2F
2773 // (Form("hTrackMatchedDEtaDPhiTrackPt%s",cutTM[i].Data()),
2774 // Form("#Delta #eta vs #Delta #varphi of cluster-track #it{p}_{T}^{track} > 0.5 GeV/#it{c}, %s",cutTM[i].Data()),
2775 // nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2776 // fhTrackMatchedDEtaDPhiTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2777 // fhTrackMatchedDEtaDPhiTrackPt[i]->SetXTitle("#Delta #eta");
2778 
2780  (Form("hTrackMatchedDEtaPosTrackPt%s",cutTM[i].Data()),
2781  Form("#Delta #eta of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2782  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2783  fhTrackMatchedDEtaPosTrackPt[i]->SetYTitle("#Delta #eta");
2784  fhTrackMatchedDEtaPosTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2785 
2787  (Form("hTrackMatchedDPhiPosTrackPt%s",cutTM[i].Data()),
2788  Form("#Delta #varphi of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2789  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2790  fhTrackMatchedDPhiPosTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2791  fhTrackMatchedDPhiPosTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2792 
2794  (Form("hTrackMatchedDEtaDPhiPosTrackPt%s",cutTM[i].Data()),
2795  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{p}_{T}^{track} > 0.5 GeV/#it{c}, %s",cutTM[i].Data()),
2796  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2797  fhTrackMatchedDEtaDPhiPosTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2798  fhTrackMatchedDEtaDPhiPosTrackPt[i]->SetXTitle("#Delta #eta");
2799 
2801  (Form("hTrackMatchedDEtaNegTrackPt%s",cutTM[i].Data()),
2802  Form("#Delta #eta of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2803  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2804  fhTrackMatchedDEtaNegTrackPt[i]->SetYTitle("#Delta #eta");
2805  fhTrackMatchedDEtaNegTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2806 
2808  (Form("hTrackMatchedDPhiNegTrackPt%s",cutTM[i].Data()),
2809  Form("#Delta #varphi of cluster-track vs #it{p}_{T}^{track}, %s",cutTM[i].Data()),
2810  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2811  fhTrackMatchedDPhiNegTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2812  fhTrackMatchedDPhiNegTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2813 
2815  (Form("hTrackMatchedDEtaDPhiNegTrackPt%s",cutTM[i].Data()),
2816  Form("#Delta #eta vs #Delta #varphi of cluster-track, #it{p}_{T}^{track} > 0.5 GeV/#it{c}, %s",cutTM[i].Data()),
2817  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2818  fhTrackMatchedDEtaDPhiNegTrackPt[i]->SetYTitle("#Delta #varphi (rad)");
2819  fhTrackMatchedDEtaDPhiNegTrackPt[i]->SetXTitle("#Delta #eta");
2820 
2821  fhdEdxTrackPt[i] = new TH2F (Form("hdEdxTrackPt%s",cutTM[i].Data()),Form("matched track <dE/dx> vs track #it{p}_{T}, %s",cutTM[i].Data()),
2822  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2823  fhdEdxTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV/#it{c})");
2824  fhdEdxTrackPt[i]->SetYTitle("<dE/dx>");
2825 
2826  fhEOverPTrackPt[i] = new TH2F (Form("hEOverPTrackPt%s",cutTM[i].Data()),Form("matched track E/p vs track #it{p}_{T}, %s",cutTM[i].Data()),
2827  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2828  fhEOverPTrackPt[i]->SetXTitle("#it{p}_{T}^{track} (GeV)^{track/#it{c}}");
2829  fhEOverPTrackPt[i]->SetYTitle("E/p");
2830 
2831 // outputContainer->Add(fhTrackMatchedDEtaTrackPt[i]) ;
2832 // outputContainer->Add(fhTrackMatchedDPhiTrackPt[i]) ;
2833 // outputContainer->Add(fhTrackMatchedDEtaDPhiTrackPt[i]) ;
2834  outputContainer->Add(fhTrackMatchedDEtaPosTrackPt[i]) ;
2835  outputContainer->Add(fhTrackMatchedDPhiPosTrackPt[i]) ;
2836  outputContainer->Add(fhTrackMatchedDEtaDPhiPosTrackPt[i]) ;
2837  outputContainer->Add(fhTrackMatchedDEtaNegTrackPt[i]) ;
2838  outputContainer->Add(fhTrackMatchedDPhiNegTrackPt[i]) ;
2839  outputContainer->Add(fhTrackMatchedDEtaDPhiNegTrackPt[i]) ;
2840  outputContainer->Add(fhdEdxTrackPt[i]);
2841  outputContainer->Add(fhEOverPTrackPt[i]);
2842  }
2843 
2845  {
2846  fhTrackMatchedDEtaTRD[i] = new TH2F
2847  (Form("hTrackMatchedDEtaTRD%s",cutTM[i].Data()),
2848  Form("#Delta #eta of cluster-track vs cluster energy, SM behind TRD, %s",cutTM[i].Data()),
2849  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2850  fhTrackMatchedDEtaTRD[i]->SetYTitle("#Delta #eta");
2851  fhTrackMatchedDEtaTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2852 
2853  fhTrackMatchedDPhiTRD[i] = new TH2F
2854  (Form("hTrackMatchedDPhiTRD%s",cutTM[i].Data()),
2855  Form("#Delta #varphi of cluster-track vs cluster energy, SM behing TRD, %s",cutTM[i].Data()),
2856  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2857  fhTrackMatchedDPhiTRD[i]->SetYTitle("#Delta #varphi (rad)");
2858  fhTrackMatchedDPhiTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2859 
2860  fhEOverPTRD[i] = new TH2F
2861  (Form("hEOverPTRD%s",cutTM[i].Data()),
2862  Form("matched track E/p vs cluster E, behind TRD, %s",cutTM[i].Data()),
2863  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2864  fhEOverPTRD[i]->SetXTitle("#it{E} (GeV)");
2865  fhEOverPTRD[i]->SetYTitle("E/p");
2866 
2867  outputContainer->Add(fhTrackMatchedDEtaTRD[i]) ;
2868  outputContainer->Add(fhTrackMatchedDPhiTRD[i]) ;
2869  outputContainer->Add(fhEOverPTRD[i]);
2870  }
2871 
2872  if(IsDataMC())
2873  {
2875  (Form("hTrackMatchedDEtaMCNoOverlap%s",cutTM[i].Data()),
2876  Form("#Delta #eta of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
2877  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2878  fhTrackMatchedDEtaMCNoOverlap[i]->SetYTitle("#Delta #eta");
2879  fhTrackMatchedDEtaMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2880 
2882  (Form("hTrackMatchedDPhiMCNoOverlap%s",cutTM[i].Data()),
2883  Form("#Delta #varphi of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
2884  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2885  fhTrackMatchedDPhiMCNoOverlap[i]->SetYTitle("#Delta #varphi (rad)");
2886  fhTrackMatchedDPhiMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2887 
2888  outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[i]) ;
2889  outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[i]) ;
2891  (Form("hTrackMatchedDEtaMCOverlap%s",cutTM[i].Data()),
2892  Form("#Delta #eta of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
2893  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2894  fhTrackMatchedDEtaMCOverlap[i]->SetYTitle("#Delta #eta");
2895  fhTrackMatchedDEtaMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2896 
2898  (Form("hTrackMatchedDPhiMCOverlap%s",cutTM[i].Data()),
2899  Form("#Delta #varphi of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
2900  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2901  fhTrackMatchedDPhiMCOverlap[i]->SetYTitle("#Delta #varphi (rad)");
2902  fhTrackMatchedDPhiMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2903 
2904  outputContainer->Add(fhTrackMatchedDEtaMCOverlap[i]) ;
2905  outputContainer->Add(fhTrackMatchedDPhiMCOverlap[i]) ;
2906 
2908  (Form("hTrackMatchedDEtaMCConversion%s",cutTM[i].Data()),
2909  Form("#Delta #eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
2910  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2911  fhTrackMatchedDEtaMCConversion[i]->SetYTitle("#Delta #eta");
2912  fhTrackMatchedDEtaMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2913 
2915  (Form("hTrackMatchedDPhiMCConversion%s",cutTM[i].Data()),
2916  Form("#Delta #varphi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
2917  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2918  fhTrackMatchedDPhiMCConversion[i]->SetYTitle("#Delta #varphi (rad)");
2919  fhTrackMatchedDPhiMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2920 
2921  outputContainer->Add(fhTrackMatchedDEtaMCConversion[i]) ;
2922  outputContainer->Add(fhTrackMatchedDPhiMCConversion[i]) ;
2923 
2925  (Form("hTrackMatchedMCParticle%s",cutTM[i].Data()),
2926  Form("Origin of particle vs energy %s",cutTM[i].Data()),
2927  nptbins,ptmin,ptmax,8,0,8);
2928  fhTrackMatchedMCParticle[i]->SetXTitle("#it{E} (GeV)");
2929  //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
2930 
2931  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(1 ,"Photon");
2932  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(2 ,"Electron");
2933  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
2934  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(4 ,"Rest");
2935  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
2936  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
2937  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
2938  fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
2939 
2940  outputContainer->Add(fhTrackMatchedMCParticle[i]);
2941  }
2942  }
2943  }
2944 
2945  if(IsPileUpAnalysisOn())
2946  {
2947  TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2948 
2949  for(Int_t i = 0 ; i < 7 ; i++)
2950  {
2951  fhPtPhotonPileUp[i] = new TH1F(Form("hPtPhotonPileUp%s",pileUpName[i].Data()),
2952  Form("Selected photon #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2953  fhPtPhotonPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2954  outputContainer->Add(fhPtPhotonPileUp[i]);
2955 
2956  fhClusterTimeDiffPhotonPileUp[i] = new TH2F(Form("hClusterTimeDiffPhotonPileUp%s",pileUpName[i].Data()),
2957  Form("Photon cluster E vs #it{t}_{max}-#it{t}_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2958  nptbins,ptmin,ptmax,400,-200,200);
2959  fhClusterTimeDiffPhotonPileUp[i]->SetXTitle("#it{E} (GeV)");
2960  fhClusterTimeDiffPhotonPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
2961  outputContainer->Add(fhClusterTimeDiffPhotonPileUp[i]);
2962  }
2963 
2964  fhTimePtPhotonNoCut = new TH2F ("hTimePtPhoton_NoCut","time of photon cluster vs pT of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2965  fhTimePtPhotonNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2966  fhTimePtPhotonNoCut->SetYTitle("#it{time} (ns)");
2967  outputContainer->Add(fhTimePtPhotonNoCut);
2968 
2969  fhTimePtPhotonSPD = new TH2F ("hTimePtPhoton_SPD","time of photon cluster vs pT of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2970  fhTimePtPhotonSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2971  fhTimePtPhotonSPD->SetYTitle("#it{time} (ns)");
2972  outputContainer->Add(fhTimePtPhotonSPD);
2973 
2974  fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,20,0,20);
2975  fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2976  fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
2977  outputContainer->Add(fhTimeNPileUpVertSPD);
2978 
2979  fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 20,0,20 );
2980  fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2981  fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
2982  outputContainer->Add(fhTimeNPileUpVertTrack);
2983 
2984  fhPtPhotonNPileUpSPDVtx = new TH2F ("hPtPhoton_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
2985  nptbins,ptmin,ptmax,20,0,20);
2986  fhPtPhotonNPileUpSPDVtx->SetYTitle("# vertex ");
2987  fhPtPhotonNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2988  outputContainer->Add(fhPtPhotonNPileUpSPDVtx);
2989 
2990  fhPtPhotonNPileUpTrkVtx = new TH2F ("hPtPhoton_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
2991  nptbins,ptmin,ptmax, 20,0,20 );
2992  fhPtPhotonNPileUpTrkVtx->SetYTitle("# vertex ");
2993  fhPtPhotonNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2994  outputContainer->Add(fhPtPhotonNPileUpTrkVtx);
2995 
2996  fhPtPhotonNPileUpSPDVtxTimeCut = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2997  nptbins,ptmin,ptmax,20,0,20);
2998  fhPtPhotonNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2999  fhPtPhotonNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3000  outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut);
3001 
3002  fhPtPhotonNPileUpTrkVtxTimeCut = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
3003  nptbins,ptmin,ptmax, 20,0,20 );
3004  fhPtPhotonNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
3005  fhPtPhotonNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3006  outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut);
3007 
3008  fhPtPhotonNPileUpSPDVtxTimeCut2 = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
3009  nptbins,ptmin,ptmax,20,0,20);
3010  fhPtPhotonNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
3011  fhPtPhotonNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3012  outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut2);
3013 
3014  fhPtPhotonNPileUpTrkVtxTimeCut2 = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
3015  nptbins,ptmin,ptmax, 20,0,20 );
3016  fhPtPhotonNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
3017  fhPtPhotonNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3018  outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut2);
3019 
3020  }
3021 
3023  {
3024  for(Int_t ieta = 0; ieta < 4; ieta++)
3025  {
3026  for(Int_t iphi = 0; iphi < 3; iphi++)
3027  {
3028 // fhLam0EMCALRegion[ieta][iphi] =
3029 // new TH2F(Form("hLam0_eta%d_phi%d",ieta,iphi),
3030 // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, region eta %d, phi %d",ieta,iphi),
3031 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3032 // fhLam0EMCALRegion[ieta][iphi]->SetYTitle("#lambda_{0}^{2}");
3033 // fhLam0EMCALRegion[ieta][iphi]->SetXTitle("#it{p}_{T} (GeV)");
3034 // outputContainer->Add(fhLam0EMCALRegion[ieta][iphi]) ;
3035 //
3036 // if(GetFirstSMCoveredByTRD() >= 0)
3037 // {
3038 // fhLam0EMCALRegionTRD[ieta][iphi] =
3039 // new TH2F(Form("hLam0TRD_eta%d_phi%d",ieta,iphi),
3040 // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, region eta %d, phi %d, SM covered by TRD",ieta,iphi),
3041 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3042 // fhLam0EMCALRegionTRD[ieta][iphi]->SetYTitle("#lambda_{0}^{2}");
3043 // fhLam0EMCALRegionTRD[ieta][iphi]->SetXTitle("#it{p}_{T} (GeV)");
3044 // outputContainer->Add(fhLam0EMCALRegionTRD[ieta][iphi]) ;
3045 // } // TRD
3046 
3047  for(Int_t ism = 0; ism < fNModules; ism++)
3048  {
3049  if(ism < fFirstModule || ism > fLastModule) continue;
3050 
3051  fhLam0EMCALRegionPerSM[ieta][iphi][ism] =
3052  new TH2F(Form("hLam0_eta%d_phi%d_sm%d",ieta,iphi,ism),
3053  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, sm %d, region eta %d, phi %d",ism,ieta,iphi),
3054  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3055  fhLam0EMCALRegionPerSM[ieta][iphi][ism]->SetYTitle("#lambda_{0}^{2}");
3056  fhLam0EMCALRegionPerSM[ieta][iphi][ism]->SetXTitle("#it{p}_{T} (GeV)");
3057  outputContainer->Add(fhLam0EMCALRegionPerSM[ieta][iphi][ism]) ;
3058 
3059  fhLam1EMCALRegionPerSM[ieta][iphi][ism] =
3060  new TH2F(Form("hLam1_eta%d_phi%d_sm%d",ieta,iphi,ism),
3061  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{1}^{2}, sm %d, region eta %d, phi %d",ism,ieta,iphi),
3062  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3063  fhLam1EMCALRegionPerSM[ieta][iphi][ism]->SetYTitle("#lambda_{1}^{2}");
3064  fhLam1EMCALRegionPerSM[ieta][iphi][ism]->SetXTitle("#it{p}_{T} (GeV)");
3065  outputContainer->Add(fhLam1EMCALRegionPerSM[ieta][iphi][ism]) ;
3066  }
3067  } // iphi
3068  } // ieta
3069 
3070  Float_t ptLimit[] = {2,3,4,5,6,8,10,12};
3071  TString l0bin [] = {"0.23<#lambda^{2}_{0}<0.26","0.3<#lambda^{2}_{0}<0.4"};
3072 
3073  for(Int_t il0 = 0; il0 < 2; il0++)
3074  {
3075  for(Int_t ipt = 0; ipt < 7; ipt++)
3076  {
3077  fhEtaPhiLam0BinPtBin[il0][ipt] = new TH2F
3078  (Form("hEtaPhiLam0Bin%d_PtBin%d",il0,ipt),
3079  Form("#eta vs #varphi in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s",
3080  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3081  netabins,etamin,etamax,nphibins,phimin,phimax);
3082  fhEtaPhiLam0BinPtBin[il0][ipt]->SetYTitle("#varphi (rad)");
3083  fhEtaPhiLam0BinPtBin[il0][ipt]->SetXTitle("#eta");
3084  outputContainer->Add(fhEtaPhiLam0BinPtBin[il0][ipt]) ;
3085 
3086  fhColRowLam0BinPtBin[il0][ipt] = new TH2F
3087  (Form("hColRowLam0Bin%d_PtBin%d",il0,ipt),
3088  Form("row vs column in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s, w > 0",
3089  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3090  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3091  fhColRowLam0BinPtBin[il0][ipt]->SetYTitle("row");
3092  fhColRowLam0BinPtBin[il0][ipt]->SetXTitle("column");
3093  outputContainer->Add(fhColRowLam0BinPtBin[il0][ipt]) ;
3094 
3095  fhColRowLam0BinPtBinWeighted[il0][ipt] = new TH2F
3096  (Form("hColRowLam0Bin%d_PtBin%dWeighted",il0,ipt),
3097  Form("cluster cell row vs column weighted in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s",
3098  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3099  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3100  fhColRowLam0BinPtBinWeighted[il0][ipt]->SetYTitle("row");
3101  fhColRowLam0BinPtBinWeighted[il0][ipt]->SetXTitle("column");
3102  outputContainer->Add(fhColRowLam0BinPtBinWeighted[il0][ipt]) ;
3103 
3104  fhColRowLam0BinPtBinLargeTime[il0][ipt] = new TH2F
3105  (Form("hColRowLam0Bin%d_PtBin%d_LargeTimeInClusterCell",il0,ipt),
3106  Form("row vs column in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s, |t| > 50 ns, w > 0",
3107  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3108  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3109  fhColRowLam0BinPtBinLargeTime[il0][ipt]->SetYTitle("row");
3110  fhColRowLam0BinPtBinLargeTime[il0][ipt]->SetXTitle("column");
3111  outputContainer->Add(fhColRowLam0BinPtBinLargeTime[il0][ipt]) ;
3112 
3113 // fhEtaPhiLam0BinPtBinSMShared[il0][ipt] = new TH2F
3114 // (Form("hEtaPhiLam0Bin%d_PtBin%d_SMShared",il0,ipt),
3115 // Form("#eta vs #varphi in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s",
3116 // ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3117 // netabins,etamin,etamax,nphibins,phimin,phimax);
3118 // fhEtaPhiLam0BinPtBinSMShared[il0][ipt]->SetYTitle("#varphi (rad)");
3119 // fhEtaPhiLam0BinPtBinSMShared[il0][ipt]->SetXTitle("#eta");
3120 // outputContainer->Add(fhEtaPhiLam0BinPtBinSMShared[il0][ipt]) ;
3121 //
3122 // fhColRowLam0BinPtBinSMShared[il0][ipt] = new TH2F
3123 // (Form("hColRowLam0Bin%d_PtBin%d_SMShared",il0,ipt),
3124 // Form("row vs column in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, %s, w > 0",
3125 // ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3126 // ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3127 // fhColRowLam0BinPtBinSMShared[il0][ipt]->SetYTitle("row");
3128 // fhColRowLam0BinPtBinSMShared[il0][ipt]->SetXTitle("column");
3129 // outputContainer->Add(fhColRowLam0BinPtBinSMShared[il0][ipt]) ;
3130 
3131  fhEtaPhiLargeTimeInClusterCell[il0][ipt] = new TH2F
3132  (Form("hEtaPhiLam0Bin%d_PtBin%d_LargeTimeInClusterCell",il0,ipt),
3133  Form("#eta vs #varphi in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, t > 50 ns, %s",
3134  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3135  netabins,etamin,etamax,nphibins,phimin,phimax);
3136  fhEtaPhiLargeTimeInClusterCell[il0][ipt]->SetYTitle("#varphi (rad)");
3137  fhEtaPhiLargeTimeInClusterCell[il0][ipt]->SetXTitle("#eta");
3138  outputContainer->Add(fhEtaPhiLargeTimeInClusterCell[il0][ipt]) ;
3139 
3140  fhCellClusterIndexEAndTime[il0][ipt] = new TH2F
3141  (Form("hCellClusterIndexEAndTimeLam0Bin%d_PtBin%d",il0,ipt),
3142  Form("#it{t}_{cell} vs cell index (E sorted) in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, %s",
3143  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3144  ntimebins,timemin,timemax,30,0,30);
3145  fhCellClusterIndexEAndTime[il0][ipt] ->SetYTitle("Index (sorted #it{E})");
3146  fhCellClusterIndexEAndTime[il0][ipt] ->SetXTitle("#it{time}_{cell} (ns)");
3147  outputContainer->Add(fhCellClusterIndexEAndTime[il0][ipt] ) ;
3148 
3149  fhCellClusterEAndTime[il0][ipt] = new TH2F
3150  (Form("hCellClusterEAndTimeLam0Bin%d_PtBin%d",il0,ipt),
3151  Form("#it{E}_{cell} vs #it{t}_{cell} in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, %s",
3152  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3153  ntimebins,timemin,timemax,100,0,5);
3154  fhCellClusterEAndTime[il0][ipt] ->SetYTitle("#it{E}_{cell} (GeV)");
3155  fhCellClusterEAndTime[il0][ipt] ->SetXTitle("#it{time}_{cell} (ns)");
3156  outputContainer->Add(fhCellClusterEAndTime[il0][ipt] ) ;
3157 
3158  fhCellClusterEFracAndTime[il0][ipt] = new TH2F
3159  (Form("hCellClusterEFracAndTimeLam0Bin%d_PtBin%d",il0,ipt),
3160  Form("#it{E}_{cell} vs #it{t}_{cell} in #it{p}_{T}=[%2.1f,%2.1f] GeV/#it{c}, w > 0, %s",
3161  ptLimit[ipt],ptLimit[ipt+1],l0bin[il0].Data()),
3162  ntimebins,timemin,timemax,100,0,1);
3163  fhCellClusterEFracAndTime[il0][ipt] ->SetYTitle("#it{E}_{cell} (GeV)");
3164  fhCellClusterEFracAndTime[il0][ipt] ->SetXTitle("#it{time}_{cell} (ns)");
3165  outputContainer->Add(fhCellClusterEFracAndTime[il0][ipt] ) ;
3166  } // pt bin
3167 
3168  for(Int_t ism = 0; ism < fNModules; ism++)
3169  {
3170  if(ism < fFirstModule || ism > fLastModule) continue;
3171 
3172  fhLam1Lam0BinPerSM[il0][ism] = new TH2F
3173  (Form("hLam1Lam0Bin%d_sm%d",il0,ism),
3174  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, %s",ism,l0bin[il0].Data()),
3175  nptbins,ptmin,ptmax,40,0,0.4);
3176  fhLam1Lam0BinPerSM[il0][ism]->SetYTitle("#lambda^{2}_{1}");
3177  fhLam1Lam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3178  outputContainer->Add(fhLam1Lam0BinPerSM[il0][ism]) ;
3179 
3180  fhTimeLam0BinPerSM[il0][ism] = new TH2F
3181  (Form("hTimeLam0Bin%d_sm%d",il0,ism),
3182  Form("#it{p}_{T} vs cluster cell time in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3183  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3184  fhTimeLam0BinPerSM[il0][ism]->SetYTitle("cell time (ns)");
3185  fhTimeLam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3186  outputContainer->Add(fhTimeLam0BinPerSM[il0][ism]) ;
3187 
3188  fhTimeLam0BinPerSMWeighted[il0][ism] = new TH2F
3189  (Form("hTimeLam0Bin%d_sm%d_Weighted",il0,ism),
3190  Form("#it{p}_{T} vs cluster cell time weighted in sm %d, %s",ism,l0bin[il0].Data()),
3191  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3192  fhTimeLam0BinPerSMWeighted[il0][ism]->SetYTitle("cell time (ns)");
3193  fhTimeLam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3194  outputContainer->Add(fhTimeLam0BinPerSMWeighted[il0][ism]) ;
3195 
3196  fhDTimeLam0BinPerSM[il0][ism] = new TH2F
3197  (Form("hDTimeLam0Bin%d_sm%d",il0,ism),
3198  Form("#it{p}_{T} vs t_{cluster}-t_{cell} in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3199  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3200  fhDTimeLam0BinPerSM[il0][ism]->SetYTitle("t_{cluster}-t_{cell} (ns)");
3201  fhDTimeLam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3202  outputContainer->Add(fhDTimeLam0BinPerSM[il0][ism]) ;
3203 
3204  fhDTimeLam0BinPerSMWeighted[il0][ism] = new TH2F
3205  (Form("hDTimeLam0Bin%d_sm%d_Weighted",il0,ism),
3206  Form("#it{p}_{T} vs t_{cluster}-t_{cell} weighted in sm %d, %s",ism,l0bin[il0].Data()),
3207  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3208  fhDTimeLam0BinPerSMWeighted[il0][ism]->SetYTitle("t_{cluster}-t_{cell} (ns)");
3209  fhDTimeLam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3210  outputContainer->Add(fhDTimeLam0BinPerSMWeighted[il0][ism]) ;
3211 
3212  fhCellClusterEFracLam0BinPerSM[il0][ism] = new TH2F
3213  (Form("hCellClusterEFracLam0Bin%d_sm%d",il0,ism),
3214  Form("#it{p}_{T} vs cell E / cluster E in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3215  nptbins,ptmin,ptmax,100,0,1);
3216  fhCellClusterEFracLam0BinPerSM[il0][ism]->SetYTitle("cell E / cluster E ");
3217  fhCellClusterEFracLam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3218  outputContainer->Add(fhCellClusterEFracLam0BinPerSM[il0][ism]) ;
3219 
3220 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism] = new TH2F
3221 // (Form("hCellClusterEFracLam0Bin%d_sm%d_Weighted",il0,ism),
3222 // Form("#it{p}_{T} vs cell E / cluster E weighted in sm %d, %s",ism,l0bin[il0].Data()),
3223 // nptbins,ptmin,ptmax,100,0,1);
3224 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism]->SetYTitle("cell E / cluster E ");
3225 // fhCellClusterEFracLam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3226 // outputContainer->Add(fhCellClusterEFracLam0BinPerSMWeighted[il0][ism]) ;
3227 
3229  (Form("hCellClusterEFracLam0Bin%d_sm%d_LargeTimeInClusterCell",il0,ism),
3230  Form("#it{p}_{T} vs cell E / cluster E in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3231  nptbins,ptmin,ptmax,100,0,1);
3232  fhCellClusterEFracLam0BinPerSMLargeTime[il0][ism]->SetYTitle("cell E / cluster E ");
3233  fhCellClusterEFracLam0BinPerSMLargeTime[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3234  outputContainer->Add(fhCellClusterEFracLam0BinPerSMLargeTime[il0][ism]) ;
3235 
3237  (Form("hCellClusterEFracLam0Bin%d_sm%d_LargeTimeInClusterCell_Total",il0,ism),
3238  Form("#it{p}_{T} vs cell E / cluster E in sm %d, %s, w > 0, all |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3239  nptbins,ptmin,ptmax,100,0,1);
3240  fhCellClusterEFracLam0BinPerSMLargeTimeTotal[il0][ism]->SetYTitle("cell E / cluster E ");
3241  fhCellClusterEFracLam0BinPerSMLargeTimeTotal[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3242  outputContainer->Add(fhCellClusterEFracLam0BinPerSMLargeTimeTotal[il0][ism]) ;
3243 
3244  fhCellClusterELam0BinPerSM[il0][ism] = new TH2F
3245  (Form("hCellClusterELam0Bin%d_sm%d",il0,ism),
3246  Form("#it{p}_{T} vs cell E in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3247  nptbins,ptmin,ptmax,500,0,10);
3248  fhCellClusterELam0BinPerSM[il0][ism]->SetYTitle("cell E (GeV)");
3249  fhCellClusterELam0BinPerSM[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3250  outputContainer->Add(fhCellClusterELam0BinPerSM[il0][ism]) ;
3251 
3253  (Form("hCellClusterELam0Bin%d_sm%d_Weighted",il0,ism),
3254  Form("#it{p}_{T} vs cell E weighted in sm %d, %s",ism,l0bin[il0].Data()),
3255  nptbins,ptmin,ptmax,500,0,10);
3256  fhCellClusterELam0BinPerSMWeighted[il0][ism]->SetYTitle("cell E (GeV)");
3257  fhCellClusterELam0BinPerSMWeighted[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3258  outputContainer->Add(fhCellClusterELam0BinPerSMWeighted[il0][ism]) ;
3259 
3261  (Form("hCellClusterELam0Bin%d_sm%d_LargeTimeInClusterCell",il0,ism),
3262  Form("#it{p}_{T} vs cell E in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3263  nptbins,ptmin,ptmax,500,0,10);
3264  fhCellClusterELam0BinPerSMLargeTime[il0][ism]->SetYTitle("cell E (GeV)");
3265  fhCellClusterELam0BinPerSMLargeTime[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3266  outputContainer->Add(fhCellClusterELam0BinPerSMLargeTime[il0][ism]) ;
3267 
3269  (Form("hCellClusterIndexELam0Bin%d_sm%d_LargeTimeInClusterCell",il0,ism),
3270  Form("#it{p}_{T} vs cell index (E sorted) in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3271  nptbins,ptmin,ptmax,30,0,30);
3272  fhCellClusterIndexELam0BinPerSMLargeTime[il0][ism]->SetYTitle("Index");
3273  fhCellClusterIndexELam0BinPerSMLargeTime[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3274  outputContainer->Add(fhCellClusterIndexELam0BinPerSMLargeTime[il0][ism]) ;
3275 
3276  fhNCellsWithLargeTimeInCluster[il0][ism] = new TH2F
3277  (Form("hNCellsWithLargeTimeInClusterLam0Bin%d_sm%d",il0,ism),
3278  Form("#it{p}_{T} vs number of cells in sm %d, %s, w > 0, |t_{cell}| > 50 ns",ism,l0bin[il0].Data()),
3279  nptbins,ptmin,ptmax,30,0,30);
3280  fhNCellsWithLargeTimeInCluster[il0][ism]->SetYTitle("#it{n}_{cells}");
3281  fhNCellsWithLargeTimeInCluster[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3282  outputContainer->Add(fhNCellsWithLargeTimeInCluster[il0][ism]) ;
3283 
3284 // if(ism < 12)
3285 // {
3286 // fhTimeLam0BinPerSMShared[il0][ism] = new TH2F
3287 // (Form("hTimeLam0Bin%d_sm%d_SMShared",il0,ism),
3288 // Form("#it{p}_{T} vs cluster cell time in sm %d, %s, w > 0",ism,l0bin[il0].Data()),
3289 // nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3290 // fhTimeLam0BinPerSMShared[il0][ism]->SetYTitle("cell time (ns)");
3291 // fhTimeLam0BinPerSMShared[il0][ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3292 // outputContainer->Add(fhTimeLam0BinPerSMShared[il0][ism]) ;
3293 // } // Run 1 SM
3294  } // sm
3295  } // l0 bin
3296 
3297  for(Int_t ism = 0; ism < fNModules; ism++)
3298  {
3299  if(ism < fFirstModule || ism > fLastModule) continue;
3300 
3301  fhLam0PerSM[ism] = new TH2F
3302  (Form("hLam0_sm%d",ism),
3303  Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d",ism),
3304  nptbins,ptmin,ptmax,40,0,0.4);
3305  fhLam0PerSM[ism]->SetYTitle("#lambda^{2}_{0}");
3306  fhLam0PerSM[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3307  outputContainer->Add(fhLam0PerSM[ism]) ;
3308 
3309  fhLam1PerSM[ism] = new TH2F
3310  (Form("hLam1_sm%d",ism),
3311  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d",ism),
3312  nptbins,ptmin,ptmax,40,0,0.4);
3313  fhLam1PerSM[ism]->SetYTitle("#lambda^{2}_{1}");
3314  fhLam1PerSM[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3315  outputContainer->Add(fhLam1PerSM[ism]) ;
3316 
3318  (Form("hLam0_sm%d_LargeTimeInClusterCell",ism),
3319  Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d,|t_{secondary cell}| > 50 ns",ism),
3320  nptbins,ptmin,ptmax,40,0,0.4);
3321  fhLam0PerSMLargeTimeInClusterCell[ism]->SetYTitle("#lambda^{2}_{0}");
3322  fhLam0PerSMLargeTimeInClusterCell[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3323  outputContainer->Add(fhLam0PerSMLargeTimeInClusterCell[ism]) ;
3324 
3326  (Form("hLam1_sm%d_LargeTimeInClusterCell",ism),
3327  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, |t_{secondary cell}| > 50 ns",ism),
3328  nptbins,ptmin,ptmax,40,0,0.4);
3329  fhLam1PerSMLargeTimeInClusterCell[ism]->SetYTitle("#lambda^{2}_{1}");
3330  fhLam1PerSMLargeTimeInClusterCell[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3331  outputContainer->Add(fhLam1PerSMLargeTimeInClusterCell[ism]) ;
3332 
3333 // fhLam0PerSMSPDPileUp[ism] = new TH2F
3334 // (Form("hLam0_sm%d_SPDPileUp",ism),
3335 // Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d, Pile-up event SPD",ism),
3336 // nptbins,ptmin,ptmax,40,0,0.4);
3337 // fhLam0PerSMSPDPileUp[ism]->SetYTitle("#lambda^{2}_{0}");
3338 // fhLam0PerSMSPDPileUp[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3339 // outputContainer->Add(fhLam0PerSMSPDPileUp[ism]) ;
3340 //
3341 // fhLam1PerSMSPDPileUp[ism] = new TH2F
3342 // (Form("hLam1_sm%d_SPDPileUp",ism),
3343 // Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, Pile-up event SPD",ism),
3344 // nptbins,ptmin,ptmax,40,0,0.4);
3345 // fhLam1PerSMSPDPileUp[ism]->SetYTitle("#lambda^{2}_{1}");
3346 // fhLam1PerSMSPDPileUp[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3347 // outputContainer->Add(fhLam1PerSMSPDPileUp[ism]) ;
3348 
3349 // if(ism < 12)
3350 // {
3351 // fhLam0PerSMShared[ism] = new TH2F
3352 // (Form("hLam0_sm%d_SMShared",ism),
3353 // Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d, SM shared",ism),
3354 // nptbins,ptmin,ptmax,40,0,0.4);
3355 // fhLam0PerSMShared[ism]->SetYTitle("#lambda^{2}_{0}");
3356 // fhLam0PerSMShared[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3357 // outputContainer->Add(fhLam0PerSMShared[ism]) ;
3358 //
3359 // fhLam1PerSMShared[ism] = new TH2F
3360 // (Form("hLam1_sm%d_SMShared",ism),
3361 // Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, SM shared",ism),
3362 // nptbins,ptmin,ptmax,40,0,0.4);
3363 // fhLam1PerSMShared[ism]->SetYTitle("#lambda^{2}_{1}");
3364 // fhLam1PerSMShared[ism]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3365 // outputContainer->Add(fhLam1PerSMShared[ism]) ;
3366 // } // run1
3367  } // sm
3368 
3369  for(Int_t ilarge = 0; ilarge < 5; ilarge++)
3370  {
3372  (Form("hLam0_NLargeTimeInClusterCell%d",ilarge),
3373  Form("#it{p}_{T} vs #lambda^{2}_{0} in sm %d,|t_{secondary cell}| > 50 ns",ilarge),
3374  nptbins,ptmin,ptmax,40,0,0.4);
3375  fhLam0PerNLargeTimeInClusterCell[ilarge]->SetYTitle("#lambda^{2}_{0}");
3376  fhLam0PerNLargeTimeInClusterCell[ilarge]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3377  outputContainer->Add(fhLam0PerNLargeTimeInClusterCell[ilarge]) ;
3378 
3380  (Form("hLam1_NLargeTimeInClusterCell%d",ilarge),
3381  Form("#it{p}_{T} vs #lambda^{2}_{1} in sm %d, |t_{secondary cell}| > 50 ns",ilarge),
3382  nptbins,ptmin,ptmax,40,0,0.4);
3383  fhLam1PerNLargeTimeInClusterCell[ilarge]->SetYTitle("#lambda^{2}_{1}");
3384  fhLam1PerNLargeTimeInClusterCell[ilarge]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3385  outputContainer->Add(fhLam1PerNLargeTimeInClusterCell[ilarge]) ;
3386  }
3387 
3388  } // regions in EMCal
3389 
3390 
3391  if(IsDataMC())
3392  {
3393  TString ptype[] = { "#gamma" , "#gamma_{#pi decay}" , "#gamma_{#eta decay}", "#gamma_{other decay}",
3394  "#pi^{0}" , "#eta" , "e^{#pm}" , "#gamma->e^{#pm}" ,
3395  "hadron?" , "Anti-N" , "Anti-P" ,
3396  "Neutron" , "Proton" , "#pi^{#pm}" ,
3397  "#gamma_{prompt}", "#gamma_{fragmentation}", "#gamma_{ISR}" , "String" } ;
3398 
3399  TString pname[] = { "Photon" , "PhotonPi0Decay" , "PhotonEtaDecay", "PhotonOtherDecay",
3400  "Pi0" , "Eta" , "Electron" , "Conversion" ,
3401  "Hadron" , "AntiNeutron" , "AntiProton" ,
3402  "Neutron" , "Proton" , "ChPion" ,
3403  "PhotonPrompt", "PhotonFragmentation", "PhotonISR" , "String" } ;
3404 
3405  for(Int_t i = 0; i < fNOriginHistograms; i++)
3406  {
3407  fhMCE[i] = new TH1F(Form("hE_MC%s",pname[i].Data()),
3408  Form("cluster from %s : E ",ptype[i].Data()),
3409  nptbins,ptmin,ptmax);
3410  fhMCE[i]->SetXTitle("#it{E} (GeV)");
3411  outputContainer->Add(fhMCE[i]) ;
3412 
3413  fhMCPt[i] = new TH1F(Form("hPt_MC%s",pname[i].Data()),
3414  Form("cluster from %s : #it{p}_{T} ",ptype[i].Data()),
3415  nptbins,ptmin,ptmax);
3416  fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3417  outputContainer->Add(fhMCPt[i]) ;
3418 
3419  fhMCEta[i] = new TH2F(Form("hEta_MC%s",pname[i].Data()),
3420  Form("cluster from %s : #eta ",ptype[i].Data()),
3421  nptbins,ptmin,ptmax,netabins,etamin,etamax);
3422  fhMCEta[i]->SetYTitle("#eta");
3423  fhMCEta[i]->SetXTitle("#it{E} (GeV)");
3424  outputContainer->Add(fhMCEta[i]) ;
3425 
3426  fhMCPhi[i] = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
3427  Form("cluster from %s : #varphi ",ptype[i].Data()),
3428  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3429  fhMCPhi[i]->SetYTitle("#varphi (rad)");
3430  fhMCPhi[i]->SetXTitle("#it{E} (GeV)");
3431  outputContainer->Add(fhMCPhi[i]) ;
3432 
3433 
3434  fhMCDeltaE[i] = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
3435  Form("MC - Reco E from %s",pname[i].Data()),
3436  nptbins,ptmin,ptmax, 200,-50,50);
3437  fhMCDeltaE[i]->SetYTitle("#Delta #it{E} (GeV)");
3438  fhMCDeltaE[i]->SetXTitle("#it{E} (GeV)");
3439  outputContainer->Add(fhMCDeltaE[i]);
3440 
3441  fhMCDeltaPt[i] = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
3442  Form("MC - Reco #it{p}_{T} from %s",pname[i].Data()),
3443  nptbins,ptmin,ptmax, 200,-50,50);
3444  fhMCDeltaPt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
3445  fhMCDeltaPt[i]->SetYTitle("#Delta #it{p}_{T} (GeV/#it{c})");
3446  outputContainer->Add(fhMCDeltaPt[i]);
3447 
3448  fhMC2E[i] = new TH2F (Form("h2E_MC%s",pname[i].Data()),
3449  Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
3450  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3451  fhMC2E[i]->SetXTitle("#it{E}_{rec} (GeV)");
3452  fhMC2E[i]->SetYTitle("#it{E}_{gen} (GeV)");
3453  outputContainer->Add(fhMC2E[i]);
3454 
3455  fhMC2Pt[i] = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
3456  Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
3457  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3458  fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
3459  fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/#it{c})");
3460  outputContainer->Add(fhMC2Pt[i]);
3461  }
3462 
3463  TString pptype[] = { "#gamma" , "#gamma_{#pi decay}" ,
3464  "#gamma_{#eta decay}", "#gamma_{other decay}" ,
3465  "#gamma_{prompt}" , "#gamma_{fragmentation}", "#gamma_{ISR}" } ;
3466 
3467  TString ppname[] = { "Photon" , "PhotonPi0Decay" ,
3468  "PhotonEtaDecay", "PhotonOtherDecay" ,
3469  "PhotonPrompt" , "PhotonFragmentation", "PhotonISR" } ;
3470 
3471  for(Int_t i = 0; i < fNPrimaryHistograms; i++)
3472  {
3473  fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
3474  Form("primary photon %s : E ",pptype[i].Data()),
3475  nptbins,ptmin,ptmax);
3476  fhEPrimMC[i]->SetXTitle("#it{E} (GeV)");
3477  outputContainer->Add(fhEPrimMC[i]) ;
3478 
3479  fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
3480  Form("primary photon %s : #it{p}_{T} ",pptype[i].Data()),
3481  nptbins,ptmin,ptmax);
3482  fhPtPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3483  outputContainer->Add(fhPtPrimMC[i]) ;
3484 
3485  fhYPrimMC[i] = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
3486  Form("primary photon %s : Rapidity ",pptype[i].Data()),
3487  nptbins,ptmin,ptmax,200,-2,2);
3488  fhYPrimMC[i]->SetYTitle("Rapidity");
3489  fhYPrimMC[i]->SetXTitle("#it{E} (GeV)");
3490  outputContainer->Add(fhYPrimMC[i]) ;
3491 
3492  fhEtaPrimMC[i] = new TH2F(Form("hEtaPrim_MC%s",ppname[i].Data()),
3493  Form("primary photon %s : #eta",pptype[i].Data()),
3494  nptbins,ptmin,ptmax,200,-2,2);
3495  fhEtaPrimMC[i]->SetYTitle("#eta");
3496  fhEtaPrimMC[i]->SetXTitle("#it{E} (GeV)");
3497  outputContainer->Add(fhEtaPrimMC[i]) ;
3498 
3499  fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
3500  Form("primary photon %s : #varphi ",pptype[i].Data()),
3501  nptbins,ptmin,ptmax,nphibins,0,TMath::TwoPi());
3502  fhPhiPrimMC[i]->SetYTitle("#varphi (rad)");
3503  fhPhiPrimMC[i]->SetXTitle("#it{E} (GeV)");
3504  outputContainer->Add(fhPhiPrimMC[i]) ;
3505 
3506 
3507  fhEPrimMCAcc[i] = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
3508  Form("primary photon %s in acceptance: E ",pptype[i].Data()),
3509  nptbins,ptmin,ptmax);
3510  fhEPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3511  outputContainer->Add(fhEPrimMCAcc[i]) ;
3512 
3513  fhPtPrimMCAcc[i] = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
3514  Form("primary photon %s in acceptance: #it{p}_{T} ",pptype[i].Data()),
3515  nptbins,ptmin,ptmax);
3516  fhPtPrimMCAcc[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3517  outputContainer->Add(fhPtPrimMCAcc[i]) ;
3518 
3519  fhYPrimMCAcc[i] = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
3520  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
3521  nptbins,ptmin,ptmax,100,-1,1);
3522  fhYPrimMCAcc[i]->SetYTitle("Rapidity");
3523  fhYPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3524  outputContainer->Add(fhYPrimMCAcc[i]) ;
3525 
3526  fhEtaPrimMCAcc[i] = new TH2F(Form("hEtaPrimAcc_MC%s",ppname[i].Data()),
3527  Form("primary photon %s in acceptance: #eta ",pptype[i].Data()),
3528  nptbins,ptmin,ptmax,netabins,etamin,etamax);
3529  fhEtaPrimMCAcc[i]->SetYTitle("#eta");
3530  fhEtaPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3531  outputContainer->Add(fhEtaPrimMCAcc[i]) ;
3532 
3533  fhPhiPrimMCAcc[i] = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
3534  Form("primary photon %s in acceptance: #varphi ",pptype[i].Data()),
3535  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3536  fhPhiPrimMCAcc[i]->SetYTitle("#varphi (rad)");
3537  fhPhiPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3538  outputContainer->Add(fhPhiPrimMCAcc[i]) ;
3539  }
3540 
3541  if(fFillSSHistograms)
3542  {
3543  TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
3544 
3545  TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
3546 
3547  for(Int_t i = 0; i < 6; i++)
3548  {
3549  fhMCELambda0[i] = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
3550  Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
3551  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3552  fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
3553  fhMCELambda0[i]->SetXTitle("#it{E} (GeV)");
3554  outputContainer->Add(fhMCELambda0[i]) ;
3555 
3556  fhMCPtLambda0[i] = new TH2F(Form("hPtLambda0_MC%s",pnamess[i].Data()),
3557  Form("cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptypess[i].Data()),
3558  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3559  fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
3560  fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3561  outputContainer->Add(fhMCPtLambda0[i]) ;
3562 
3563  if(!GetReader()->IsEmbeddedClusterSelectionOn())
3564  {
3565  for(Int_t iover = 0; iover < 3; iover++)
3566  {
3567  fhMCPtLambda0Overlaps[i][iover] = new TH2F(Form("hPtLambda0_MC%s_Overlap%d",pnamess[i].Data(),iover),
3568  Form("cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, N Overlaps = %d",ptypess[i].Data(),iover),
3569  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3570  fhMCPtLambda0Overlaps[i][iover]->SetYTitle("#lambda_{0}^{2}");
3571  fhMCPtLambda0Overlaps[i][iover]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3572  outputContainer->Add(fhMCPtLambda0Overlaps[i][iover]) ;
3573  }
3574  }
3575 
3576  fhMCELambda1[i] = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
3577  Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
3578  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3579  fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
3580  fhMCELambda1[i]->SetXTitle("#it{E} (GeV)");
3581  outputContainer->Add(fhMCELambda1[i]) ;
3582 
3583  fhMCNCellsE[i] = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
3584  Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
3585  nptbins,ptmin,ptmax, nbins,nmin,nmax);
3586  fhMCNCellsE[i]->SetXTitle("#it{E} (GeV)");
3587  fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
3588  outputContainer->Add(fhMCNCellsE[i]);
3589 
3591  {
3592  fhMCEDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
3593  Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
3594  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3595  fhMCEDispersion[i]->SetYTitle("D^{2}");
3596  fhMCEDispersion[i]->SetXTitle("#it{E} (GeV)");
3597  outputContainer->Add(fhMCEDispersion[i]) ;
3598 
3599  fhMCMaxCellDiffClusterE[i] = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
3600  Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
3601  nptbins,ptmin,ptmax, 500,0,1.);
3602  fhMCMaxCellDiffClusterE[i]->SetXTitle("#it{E}_{cluster} (GeV) ");
3603  fhMCMaxCellDiffClusterE[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3604  outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
3605 
3606  fhMCLambda0vsClusterMaxCellDiffE0[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
3607  Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
3608  ssbins,ssmin,ssmax,500,0,1.);
3609  fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
3610  fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3611  outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
3612 
3613  fhMCLambda0vsClusterMaxCellDiffE2[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
3614  Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
3615  ssbins,ssmin,ssmax,500,0,1.);
3616  fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
3617  fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3618  outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
3619 
3620  fhMCLambda0vsClusterMaxCellDiffE6[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
3621  Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
3622  ssbins,ssmin,ssmax,500,0,1.);
3623  fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
3624  fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3625  outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
3626 
3627  fhMCNCellsvsClusterMaxCellDiffE0[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
3628  Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
3629  nbins/5,nmin,nmax/5,500,0,1.);
3630  fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
3631  fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3632  outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
3633 
3634  fhMCNCellsvsClusterMaxCellDiffE2[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
3635  Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
3636  nbins/5,nmin,nmax/5,500,0,1.);
3637  fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
3638  fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3639  outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
3640 
3641  fhMCNCellsvsClusterMaxCellDiffE6[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
3642  Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
3643  nbins/5,nmin,nmax/5,500,0,1.);
3644  fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
3645  fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("#it{E} (GeV)");
3646  outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
3647 
3648  if(GetCalorimeter()==kEMCAL)
3649  {
3650  fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
3651  Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
3652  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3653  fhMCEDispEta[i]->SetXTitle("#it{E} (GeV)");
3654  fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
3655  outputContainer->Add(fhMCEDispEta[i]);
3656 
3657  fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
3658  Form("cluster from %s : #sigma^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i} - <#varphi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
3659  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3660  fhMCEDispPhi[i]->SetXTitle("#it{E} (GeV)");
3661  fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3662  outputContainer->Add(fhMCEDispPhi[i]);
3663 
3664  fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
3665  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()),
3666  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
3667  fhMCESumEtaPhi[i]->SetXTitle("#it{E} (GeV)");
3668  fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #varphi}");
3669  outputContainer->Add(fhMCESumEtaPhi[i]);
3670 
3671  fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
3672  Form("cluster from %s : #sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),
3673  nptbins,ptmin,ptmax,200,-10,10);
3674  fhMCEDispEtaPhiDiff[i]->SetXTitle("#it{E} (GeV)");
3675  fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#varphi #varphi}-#sigma^{2}_{#eta #eta}");
3676  outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
3677 
3678  fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
3679  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()),
3680  nptbins,ptmin,ptmax, 200,-1,1);
3681  fhMCESphericity[i]->SetXTitle("#it{E} (GeV)");
3682  fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi})");
3683  outputContainer->Add(fhMCESphericity[i]);
3684 
3685  for(Int_t ie = 0; ie < 7; ie++)
3686  {
3687  fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
3688  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]),
3689  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3690  fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
3691  fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3692  outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
3693 
3694  fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pnamess[i].Data()),
3695  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]),
3696  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3697  fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
3698  fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3699  outputContainer->Add(fhMCLambda0DispEta[ie][i]);
3700 
3701  fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
3702  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]),
3703  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3704  fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
3705  fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
3706  outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
3707  }
3708  }
3709  }
3710  }// loop
3711 
3712 // if(!GetReader()->IsEmbeddedClusterSelectionOn())
3713 // {
3714 // fhMCPhotonELambda0NoOverlap = new TH2F("hELambda0_MCPhoton_NoOverlap",
3715 // "cluster from Photon : E vs #lambda_{0}^{2}",
3716 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3717 // fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
3718 // fhMCPhotonELambda0NoOverlap->SetXTitle("#it{E} (GeV)");
3719 // outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
3720 //
3721 // fhMCPhotonELambda0TwoOverlap = new TH2F("hELambda0_MCPhoton_TwoOverlap",
3722 // "cluster from Photon : E vs #lambda_{0}^{2}",
3723 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3724 // fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
3725 // fhMCPhotonELambda0TwoOverlap->SetXTitle("#it{E} (GeV)");
3726 // outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
3727 //
3728 // fhMCPhotonELambda0NOverlap = new TH2F("hELambda0_MCPhoton_NOverlap",
3729 // "cluster from Photon : E vs #lambda_{0}^{2}",
3730 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3731 // fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
3732 // fhMCPhotonELambda0NOverlap->SetXTitle("#it{E} (GeV)");
3733 // outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
3734 // } // No embedding
3735 
3736  if(GetReader()->IsEmbeddedClusterSelectionOn())
3737  {
3738  fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
3739  "Energy Fraction of embedded signal versus cluster energy",
3740  nptbins,ptmin,ptmax,100,0.,1.);
3741  fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
3742  fhEmbeddedSignalFractionEnergy->SetXTitle("#it{E} (GeV)");
3743  outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
3744 
3745  fhEmbedPhotonELambda0FullSignal = new TH2F("hELambda0_EmbedPhoton_FullSignal",
3746  "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
3747  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3748  fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
3749  fhEmbedPhotonELambda0FullSignal->SetXTitle("#it{E} (GeV)");
3750  outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
3751 
3752  fhEmbedPhotonELambda0MostlySignal = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
3753  "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
3754  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3755  fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
3756  fhEmbedPhotonELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
3757  outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
3758 
3759  fhEmbedPhotonELambda0MostlyBkg = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
3760  "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
3761  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3762  fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
3763  fhEmbedPhotonELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
3764  outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
3765 
3766  fhEmbedPhotonELambda0FullBkg = new TH2F("hELambda0_EmbedPhoton_FullBkg",
3767  "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
3768  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3769  fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
3770  fhEmbedPhotonELambda0FullBkg->SetXTitle("#it{E} (GeV)");
3771  outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
3772 
3773  fhEmbedPi0ELambda0FullSignal = new TH2F("hELambda0_EmbedPi0_FullSignal",
3774  "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
3775  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3776  fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
3777  fhEmbedPi0ELambda0FullSignal->SetXTitle("#it{E} (GeV)");
3778  outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
3779 
3780  fhEmbedPi0ELambda0MostlySignal = new TH2F("hELambda0_EmbedPi0_MostlySignal",
3781  "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
3782  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3783  fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
3784  fhEmbedPi0ELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
3785  outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
3786 
3787  fhEmbedPi0ELambda0MostlyBkg = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
3788  "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
3789  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3790  fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
3791  fhEmbedPi0ELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
3792  outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
3793 
3794  fhEmbedPi0ELambda0FullBkg = new TH2F("hELambda0_EmbedPi0_FullBkg",
3795  "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
3796  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3797  fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
3798  fhEmbedPi0ELambda0FullBkg->SetXTitle("#it{E} (GeV)");
3799  outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
3800  }// embedded histograms
3801  }// Fill SS MC histograms
3802 
3803 
3805  {
3806  fhMCConversionVertex = new TH2F("hMCPhotonConversionVertex","cluster from converted photon, #it{p}_{T} vs vertex distance",
3807  nptbins,ptmin,ptmax,500,0,500);
3808  fhMCConversionVertex->SetYTitle("#it{R} (cm)");
3809  fhMCConversionVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3810  outputContainer->Add(fhMCConversionVertex) ;
3811 
3812  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
3813  {
3814  fhMCConversionVertexTRD = new TH2F("hMCPhotonConversionVertexTRD","cluster from converted photon, #it{p}_{T} vs vertex distance, SM covered by TRD",
3815  nptbins,ptmin,ptmax,500,0,500);
3816  fhMCConversionVertexTRD->SetYTitle("#it{R} (cm)");
3817  fhMCConversionVertexTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3818  outputContainer->Add(fhMCConversionVertexTRD) ;
3819  }
3820 
3821  if(fFillSSHistograms)
3822  {
3823  TString region[] = {"ITS","TPC","TRD","TOF","Top EMCal","In EMCal"};
3824  for(Int_t iR = 0; iR < 6; iR++)
3825  {
3826  fhMCConversionLambda0Rcut[iR] = new TH2F(Form("hMCPhotonConversionLambda0_R%d",iR),
3827  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s",region[iR].Data()),
3828  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3829  fhMCConversionLambda0Rcut[iR]->SetYTitle("#lambda_{0}^{2}");
3830  fhMCConversionLambda0Rcut[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3831  outputContainer->Add(fhMCConversionLambda0Rcut[iR]) ;
3832 
3833  fhMCConversionLambda1Rcut[iR] = new TH2F(Form("hMCPhotonConversionLambda1_R%d",iR),
3834  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{1}^{2}, conversion in %s",region[iR].Data()),
3835  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3836  fhMCConversionLambda1Rcut[iR]->SetYTitle("#lambda_{1}^{2}");
3837  fhMCConversionLambda1Rcut[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3838  outputContainer->Add(fhMCConversionLambda1Rcut[iR]) ;
3839  } // R cut
3840 
3841 
3842  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
3843  {
3844  for(Int_t iR = 0; iR < 6; iR++)
3845  {
3846  fhMCConversionLambda0RcutTRD[iR] = new TH2F(Form("hMCPhotonConversionLambda0TRD_R%d",iR),
3847  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s, SM covered by TRD",region[iR].Data()),
3848  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3849  fhMCConversionLambda0RcutTRD[iR]->SetYTitle("#lambda_{0}^{2}");
3850  fhMCConversionLambda0RcutTRD[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3851  outputContainer->Add(fhMCConversionLambda0RcutTRD[iR]) ;
3852 
3853  fhMCConversionLambda1RcutTRD[iR] = new TH2F(Form("hMCPhotonConversionLambda1TRD_R%d",iR),
3854  Form("cluster from converted photon, #it{p}_{T} vs #lambda_{1}^{2}, conversion in %s, SM covered by TRD",region[iR].Data()),
3855  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3856  fhMCConversionLambda1RcutTRD[iR]->SetYTitle("#lambda_{1}^{2}");
3857  fhMCConversionLambda1RcutTRD[iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3858  outputContainer->Add(fhMCConversionLambda1RcutTRD[iR]) ;
3859  } // R cut
3860  }
3861 
3862  // if(GetCalorimeter() == kEMCAL && fFillEMCALRegionSSHistograms)
3863  // {
3864  // for(Int_t ieta = 0; ieta < 4; ieta++)
3865  // {
3866  // for(Int_t iphi = 0; iphi < 3; iphi++)
3867  // {
3868  // for(Int_t iR = 0; iR < 6; iR++)
3869  // {
3870  // fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR] =
3871  // new TH2F(Form("hMCPhotonConversionLambda0_R%d_eta%d_phi%d",iR,ieta,iphi),
3872  // 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),
3873  // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3874  // fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR]->SetYTitle("#lambda_{0}^{2}");
3875  // fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3876  // outputContainer->Add(fhLam0EMCALRegionMCConvRcut[ieta][iphi][iR]) ;
3877  //
3878  // if(GetFirstSMCoveredByTRD() >= 0)
3879  // {
3880  // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR] =
3881  // new TH2F(Form("hMCPhotonConversionLambda0TRD_R%d_eta%d_phi%d",iR,ieta,iphi),
3882  // Form("cluster from converted photon, #it{p}_{T} vs #lambda_{0}^{2}, conversion in %s, region eta %d, phi %d, SM covered by TRD",region[iR].Data(),ieta,iphi),
3883  // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3884  // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR]->SetYTitle("#lambda_{0}^{2}");
3885  // fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3886  // outputContainer->Add(fhLam0EMCALRegionTRDMCConvRcut[ieta][iphi][iR]) ;
3887  // } // TRD
3888  //
3889  // } // iR
3890  // } // iphi
3891  // } // ieta
3892  // } // regions in EMCal
3893  } // shower shape
3894  } // conversion vertex
3895  } // Histos with MC
3896 
3898  {
3899  for(Int_t ie=0; ie<fNEBinCuts; ie++)
3900  {
3901  fhEBinClusterEtaPhi[ie] = new TH2F
3902  (Form("hEBin%d_Cluster_EtaPhi",ie),
3903  Form("#eta vs #varphi, cluster, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
3904  netabins,etamin,etamax,nphibins,phimin,phimax);
3905  fhEBinClusterEtaPhi[ie]->SetYTitle("#varphi (rad)");
3906  fhEBinClusterEtaPhi[ie]->SetXTitle("#eta");
3907  outputContainer->Add(fhEBinClusterEtaPhi[ie]) ;
3908 
3909  fhEBinClusterColRow[ie] = new TH2F
3910  (Form("hEBin%d_Cluster_ColRow",ie),
3911  Form("column vs row, cluster max E cell, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
3912  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3913  fhEBinClusterColRow[ie]->SetYTitle("row");
3914  fhEBinClusterColRow[ie]->SetXTitle("column");
3915  outputContainer->Add(fhEBinClusterColRow[ie]) ;
3916 
3917  fhEBinClusterEtaPhiPID[ie] = new TH2F
3918  (Form("hEBin%d_Cluster_EtaPhi_PID",ie),
3919  Form("#eta vs #varphi, cluster, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}, PID cut",fEBinCuts[ie],fEBinCuts[ie+1]),
3920  netabins,etamin,etamax,nphibins,phimin,phimax);
3921  fhEBinClusterEtaPhiPID[ie]->SetYTitle("#varphi (rad)");
3922  fhEBinClusterEtaPhiPID[ie]->SetXTitle("#eta");
3923  outputContainer->Add(fhEBinClusterEtaPhiPID[ie]) ;
3924 
3925  fhEBinClusterColRowPID[ie] = new TH2F
3926  (Form("hEBin%d_Cluster_ColRow_PID",ie),
3927  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]),
3928  ncolcell,colcellmin,colcellmax,nrowcell,rowcellmin,rowcellmax);
3929  fhEBinClusterColRowPID[ie]->SetYTitle("row");
3930  fhEBinClusterColRowPID[ie]->SetXTitle("column");
3931  outputContainer->Add(fhEBinClusterColRowPID[ie]) ;
3932  }
3933  }
3934 
3936  {
3937  TString caseTitle[] = {"","CleanCluster","MergedClusterHijingBkg","MergedClusterNotHijingBkg","MergedClusterHijingAndOtherBkg","MergedCluster"};
3938  Int_t ncases = 1;
3939  if ( IsDataMC() && IsStudyClusterOverlapsPerGeneratorOn() ) ncases = 6;
3940 
3941  for(Int_t icase = 0; icase < ncases; icase++)
3942  {
3943  fhLocalRegionClusterEtaPhi[icase] = new TH2F
3944  (Form("hLocalRegionClusterEtaPhi%s",caseTitle[icase].Data()),
3945  "cluster,#it{E} > 0.5 GeV, #eta vs #varphi",netabins,etamin,etamax,nphibins,phimin,phimax);
3946  fhLocalRegionClusterEtaPhi[icase]->SetYTitle("#varphi (rad)");
3947  fhLocalRegionClusterEtaPhi[icase]->SetXTitle("#eta");
3948  outputContainer->Add(fhLocalRegionClusterEtaPhi[icase]) ;
3949 
3950  fhLocalRegionClusterEnergySum[icase] = new TH2F
3951  (Form("hLocalRegionClusterEnergySum%s",caseTitle[icase].Data()),
3952  "Sum of cluster energy around trigger cluster #it{E} with R=0.2",
3953  nptbins,ptmin,ptmax, 200,0,100);
3954  fhLocalRegionClusterEnergySum[icase]->SetXTitle("#it{E} (GeV)");
3955  fhLocalRegionClusterEnergySum[icase]->SetYTitle("#Sigma #it{E} (GeV)");
3956  outputContainer->Add(fhLocalRegionClusterEnergySum[icase]);
3957 
3959  (Form("hLocalRegionClusterMultiplicity%s",caseTitle[icase].Data()),
3960  "Cluster multiplicity around trigger cluster #it{E} with R=0.2",
3961  nptbins,ptmin,ptmax, 200,0,200);
3962  fhLocalRegionClusterMultiplicity[icase]->SetXTitle("#it{E} (GeV)");
3963  fhLocalRegionClusterMultiplicity[icase]->SetYTitle("Multiplicity");
3964  outputContainer->Add(fhLocalRegionClusterMultiplicity[icase]);
3965 
3967  {
3969  (Form("hLocalRegionClusterEnergySumPerCentrality%s",caseTitle[icase].Data()),
3970  "Sum of cluster energy around trigger cluster vs centrality with R=0.2",
3971  100,0,100, 200,0,100);
3972  fhLocalRegionClusterEnergySumPerCentrality[icase]->SetXTitle("Centrality");
3973  fhLocalRegionClusterEnergySumPerCentrality[icase]->SetYTitle("#Sigma #it{E} (GeV)");
3974  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentrality[icase]);
3975 
3977  (Form("hLocalRegionClusterMultiplicityPerCentrality%s",caseTitle[icase].Data()),
3978  "Cluster multiplicity around trigger cluster vs centrality with R=0.2",
3979  100,0,100, 200,0,200);
3980  fhLocalRegionClusterMultiplicityPerCentrality[icase]->SetXTitle("Centrality");
3981  fhLocalRegionClusterMultiplicityPerCentrality[icase]->SetYTitle("Multiplicity");
3982  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentrality[icase]);
3983  }
3984 
3985  //
3986  // Select only single pi0 decay clusters from non hijing generator
3987  //
3988  if(IsDataMC())
3989  {
3991  (Form("hLocalRegionClusterEnergySum%s_MCPi0Decay",caseTitle[icase].Data()),
3992  "Sum of cluster energy around trigger cluster #it{E} with R=0.2",
3993  nptbins,ptmin,ptmax, 200,0,100);
3994  fhLocalRegionClusterEnergySumMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
3995  fhLocalRegionClusterEnergySumMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
3996  outputContainer->Add(fhLocalRegionClusterEnergySumMCPi0Decay[icase]);
3997 
3999  (Form("hLocalRegionClusterMultiplicity%s_MCPi0Decay",caseTitle[icase].Data()),
4000  "Cluster multiplicity around trigger cluster #it{E} with R=0.2",
4001  nptbins,ptmin,ptmax, 200,0,200);
4002  fhLocalRegionClusterMultiplicityMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4003  fhLocalRegionClusterMultiplicityMCPi0Decay[icase]->SetYTitle("Multiplicity");
4004  outputContainer->Add(fhLocalRegionClusterMultiplicityMCPi0Decay[icase]);
4005 
4007  {
4009  (Form("hLocalRegionClusterEnergySumPerCentrality%s_MCPi0Decay",caseTitle[icase].Data()),
4010  "Sum of cluster energy around trigger cluster vs centrality with R=0.2",
4011  100,0,100, 200,0,100);
4012  fhLocalRegionClusterEnergySumPerCentralityMCPi0Decay[icase]->SetXTitle("Centrality");
4013  fhLocalRegionClusterEnergySumPerCentralityMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4014  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityMCPi0Decay[icase]);
4015 
4017  (Form("hLocalRegionClusterMultiplicityPerCentrality%s_MCPi0Decay",caseTitle[icase].Data()),
4018  "Cluster multiplicity around trigger cluster vs centrality with R=0.2",
4019  100,0,100, 200,0,200);
4020  fhLocalRegionClusterMultiplicityPerCentralityMCPi0Decay[icase]->SetXTitle("Centrality");
4021  fhLocalRegionClusterMultiplicityPerCentralityMCPi0Decay[icase]->SetYTitle("Multiplicity");
4022  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentralityMCPi0Decay[icase]);
4023  }
4024  }
4025 
4027  {
4029  (Form("hLocalRegionClusterEnergySumHijing%s",caseTitle[icase].Data()),
4030  "Sum of cluster energy (HIJING) around trigger cluster #it{E} with R=0.2",
4031  nptbins,ptmin,ptmax, 200,0,100);
4032  fhLocalRegionClusterEnergySumHijing[icase]->SetXTitle("#it{E} (GeV)");
4033  fhLocalRegionClusterEnergySumHijing[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4034  outputContainer->Add(fhLocalRegionClusterEnergySumHijing[icase]);
4035 
4037  (Form("hLocalRegionClusterMultiplicityHijing%s",caseTitle[icase].Data()),
4038  "Cluster multiplicity (HIJING) around trigger cluster #it{E} with R=0.2",
4039  nptbins,ptmin,ptmax, 200,0,200);
4040  fhLocalRegionClusterMultiplicityHijing[icase]->SetXTitle("#it{E} (GeV)");
4041  fhLocalRegionClusterMultiplicityHijing[icase]->SetYTitle("Multiplicity");
4042  outputContainer->Add(fhLocalRegionClusterMultiplicityHijing[icase]);
4043 
4045  (Form("hLocalRegionClusterEnergySumAdded%s",caseTitle[icase].Data()),
4046  "Sum of cluster energy (not HIJING) around trigger cluster #it{E} with R=0.2",
4047  nptbins,ptmin,ptmax, 200,0,100);
4048  fhLocalRegionClusterEnergySumAdded[icase]->SetXTitle("#it{E} (GeV)");
4049  fhLocalRegionClusterEnergySumAdded[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4050  outputContainer->Add(fhLocalRegionClusterEnergySumAdded[icase]);
4051 
4053  (Form("hLocalRegionClusterMultiplicityAdded%s",caseTitle[icase].Data()),
4054  "Cluster multiplicity (not HIJING) around trigger cluster #it{E} with R=0.2",
4055  nptbins,ptmin,ptmax, 200,0,200);
4056  fhLocalRegionClusterMultiplicityAdded[icase]->SetXTitle("#it{E} (GeV)");
4057  fhLocalRegionClusterMultiplicityAdded[icase]->SetYTitle("Multiplicity");
4058  outputContainer->Add(fhLocalRegionClusterMultiplicityAdded[icase]);
4059 
4060 
4062  {
4064  (Form("hLocalRegionClusterEnergySumPerCentralityHijing%s",caseTitle[icase].Data()),
4065  "Sum of cluster energy (HIJING) around trigger cluster vs centrality with R=0.2",
4066  100,0,100, 200,0,100);
4067  fhLocalRegionClusterEnergySumPerCentralityHijing[icase]->SetXTitle("Centrality");
4068  fhLocalRegionClusterEnergySumPerCentralityHijing[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4069  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityHijing[icase]);
4070 
4072  (Form("hLocalRegionClusterMultiplicityPerCentralityHijing%s",caseTitle[icase].Data()),
4073  "Cluster multiplicity (HIJING) around trigger cluster vs centrality with R=0.2",
4074  100,0,100, 200,0,200);
4075  fhLocalRegionClusterMultiplicityPerCentralityHijing[icase]->SetXTitle("Centrality");
4076  fhLocalRegionClusterMultiplicityPerCentralityHijing[icase]->SetYTitle("Multiplicity");
4077  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentralityHijing[icase]);
4078 
4079 
4081  (Form("hLocalRegionClusterEnergySumPerCentralityAdded%s",caseTitle[icase].Data()),
4082  "Sum of cluster energy (not HIJING) around trigger cluster vs centrality with R=0.2",
4083  100,0,100, 200,0,100);
4084  fhLocalRegionClusterEnergySumPerCentralityAdded[icase]->SetXTitle("Centrality");
4085  fhLocalRegionClusterEnergySumPerCentralityAdded[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4086  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityAdded[icase]);
4087 
4089  (Form("hLocalRegionClusterMultiplicityPerCentralityAdded%s",caseTitle[icase].Data()),
4090  "Cluster multiplicity (not HIJING) around trigger cluster vs centrality with R=0.2",
4091  100,0,100, 200,0,200);
4092  fhLocalRegionClusterMultiplicityPerCentralityAdded[icase]->SetXTitle("Centrality");
4093  fhLocalRegionClusterMultiplicityPerCentralityAdded[icase]->SetYTitle("Multiplicity");
4094  outputContainer->Add(fhLocalRegionClusterMultiplicityPerCentralityAdded[icase]);
4095  }
4096 
4097  //
4098  // Select only single pi0 decay clusters from non hijing generator
4099  //
4100 
4102  (Form("hLocalRegionClusterEnergySumHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4103  "Sum of cluster energy (HIJING) around trigger cluster #it{E} with R=0.2",
4104  nptbins,ptmin,ptmax, 200,0,100);
4105  fhLocalRegionClusterEnergySumHijingMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4106  fhLocalRegionClusterEnergySumHijingMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4107  outputContainer->Add(fhLocalRegionClusterEnergySumHijingMCPi0Decay[icase]);
4108 
4110  (Form("hLocalRegionClusterMultiplicityHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4111  "Cluster multiplicity (HIJING) around trigger cluster #it{E} with R=0.2",
4112  nptbins,ptmin,ptmax, 200,0,200);
4113  fhLocalRegionClusterMultiplicityHijingMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4114  fhLocalRegionClusterMultiplicityHijingMCPi0Decay[icase]->SetYTitle("Multiplicity");
4115  outputContainer->Add(fhLocalRegionClusterMultiplicityHijingMCPi0Decay[icase]);
4116 
4118  (Form("hLocalRegionClusterEnergySumAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4119  "Sum of cluster energy (not HIJING) around trigger cluster #it{E} with R=0.2",
4120  nptbins,ptmin,ptmax, 200,0,100);
4121  fhLocalRegionClusterEnergySumAddedMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4122  fhLocalRegionClusterEnergySumAddedMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4123  outputContainer->Add(fhLocalRegionClusterEnergySumAddedMCPi0Decay[icase]);
4124 
4126  (Form("hLocalRegionClusterMultiplicityAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4127  "Cluster multiplicity (not HIJING) around trigger cluster #it{E} with R=0.2",
4128  nptbins,ptmin,ptmax, 200,0,200);
4129  fhLocalRegionClusterMultiplicityAddedMCPi0Decay[icase]->SetXTitle("#it{E} (GeV)");
4130  fhLocalRegionClusterMultiplicityAddedMCPi0Decay[icase]->SetYTitle("Multiplicity");
4131  outputContainer->Add(fhLocalRegionClusterMultiplicityAddedMCPi0Decay[icase]);
4132 
4133 
4135  {
4137  (Form("hLocalRegionClusterEnergySumPerCentralityHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4138  "Sum of cluster energy (HIJING) around trigger cluster vs centrality with R=0.2",
4139  100,0,100, 200,0,100);
4140  fhLocalRegionClusterEnergySumPerCentralityHijingMCPi0Decay[icase]->SetXTitle("Centrality");
4141  fhLocalRegionClusterEnergySumPerCentralityHijingMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4142  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityHijingMCPi0Decay[icase]);
4143 
4145  (Form("hLocalRegionClusterMultiplicityPerCentralityHijing%s_MCPi0Decay",caseTitle[icase].Data()),
4146  "Cluster multiplicity (HIJING) around trigger cluster vs centrality with R=0.2",
4147  100,0,100, 200,0,200);
4148  fhLocalRegionClusterMultiplicityPerCentralityHijingMCPi0Decay[icase]->SetXTitle("Centrality");
4149  fhLocalRegionClusterMultiplicityPerCentralityHijingMCPi0Decay[icase]->SetYTitle("Multiplicity");
4151 
4152 
4154  (Form("hLocalRegionClusterEnergySumPerCentralityAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4155  "Sum of cluster energy (not HIJING) around trigger cluster vs centrality with R=0.2",
4156  100,0,100, 200,0,100);
4157  fhLocalRegionClusterEnergySumPerCentralityAddedMCPi0Decay[icase]->SetXTitle("Centrality");
4158  fhLocalRegionClusterEnergySumPerCentralityAddedMCPi0Decay[icase]->SetYTitle("#Sigma #it{E} (GeV)");
4159  outputContainer->Add(fhLocalRegionClusterEnergySumPerCentralityAddedMCPi0Decay[icase]);
4160 
4162  (Form("hLocalRegionClusterMultiplicityPerCentralityAdded%s_MCPi0Decay",caseTitle[icase].Data()),
4163  "Cluster multiplicity (not HIJING) around trigger cluster vs centrality with R=0.2",
4164  100,0,100, 200,0,200);
4165  fhLocalRegionClusterMultiplicityPerCentralityAddedMCPi0Decay[icase]->SetXTitle("Centrality");
4166  fhLocalRegionClusterMultiplicityPerCentralityAddedMCPi0Decay[icase]->SetYTitle("Multiplicity");
4168  }
4169 
4170  }
4171  }
4172 
4173  if(IsDataMC())
4174  {
4176  ("hLocalRegionClusterEnergySumHijing2",
4177  "Sum of cluster energy (HIJING) around trigger cluster #it{E} with R=0.2",
4178  nptbins,ptmin,ptmax, 200,0,100);
4179  fhLocalRegionClusterEnergySumHijing2->SetXTitle("#it{E} (GeV)");
4180  fhLocalRegionClusterEnergySumHijing2->SetYTitle("#Sigma #it{E} (GeV)");
4181  outputContainer->Add(fhLocalRegionClusterEnergySumHijing2);
4182 
4184  ("hLocalRegionClusterMultiplicityHijing2",
4185  "Cluster multiplicity (HIJING) around trigger cluster #it{E} with R=0.2",
4186  nptbins,ptmin,ptmax, 200,0,200);
4187  fhLocalRegionClusterMultiplicityHijing2->SetXTitle("#it{E} (GeV)");
4188  fhLocalRegionClusterMultiplicityHijing2->SetYTitle("Multiplicity");
4189  outputContainer->Add(fhLocalRegionClusterMultiplicityHijing2);
4190 
4192  {
4194  ("hLocalRegionClusterEnergySumPerCentralityHijing2",
4195  "Sum of cluster energy (HIJING) around trigger cluster vs centrality with R=0.2",
4196  100,0,100, 200,0,100);
4197  fhLocalRegionClusterEnergySumPerCentralityHijing2->SetXTitle("Centrality");
4198  fhLocalRegionClusterEnergySumPerCentralityHijing2->SetYTitle("#Sigma #it{E} (GeV)");
4200 
4202  ("hLocalRegionClusterMultiplicityPerCentralityHijing2",
4203  "Cluster multiplicity (HIJING) around trigger cluster vs centrality with R=0.2",
4204  100,0,100, 200,0,200);
4206  fhLocalRegionClusterMultiplicityPerCentralityHijing2->SetYTitle("Multiplicity");
4208  }
4209  }
4210  // fhDistanceAddedPhotonAddedPrimarySignal = new TH2F
4211  // ("hDistanceAddedPhotonAddedPrimarySignal", "Distance added #gamma to primary particle from added generator"
4212  // ,nptbins,ptmin,ptmax,100,0,0.4);
4213  // fhDistanceAddedPhotonAddedPrimarySignal->SetYTitle("#it{R}");
4214  // fhDistanceAddedPhotonAddedPrimarySignal->SetXTitle("#it{p}_{T}");
4215  // outputContainer->Add(fhDistanceAddedPhotonAddedPrimarySignal) ;
4216  //
4217  // fhDistanceHijingPhotonAddedPrimarySignal = new TH2F
4218  // ("hDistanceHijingPhotonAddedPrimarySignal", "Distance Hijing #gamma to primary particle from added generator"
4219  // ,nptbins,ptmin,ptmax,100,0,0.4);
4220  // fhDistanceHijingPhotonAddedPrimarySignal->SetYTitle("#it{R}");
4221  // fhDistanceHijingPhotonAddedPrimarySignal->SetXTitle("#it{p}_{T}");
4222  // outputContainer->Add(fhDistanceHijingPhotonAddedPrimarySignal) ;
4223  //
4224  // fhDistanceAddedPhotonAddedSecondarySignal = new TH2F
4225  // ("hDistanceAddedPhotonAddedSecondarySignal", "Distance added #gamma to secondary particle from added generator"
4226  // ,nptbins,ptmin,ptmax,100,0,0.4);
4227  // fhDistanceAddedPhotonAddedSecondarySignal->SetYTitle("#it{R}");
4228  // fhDistanceAddedPhotonAddedSecondarySignal->SetXTitle("#it{p}_{T}");
4229  // outputContainer->Add(fhDistanceAddedPhotonAddedSecondarySignal) ;
4230  //
4231  // fhDistanceHijingPhotonAddedSecondarySignal = new TH2F
4232  // ("hDistanceHijingPhotonAddedSecondarySignal", "Distance Hijing #gamma to secondary particle from added generator"
4233  // ,nptbins,ptmin,ptmax,100,0,0.4);
4234  // fhDistanceHijingPhotonAddedSecondarySignal->SetYTitle("#it{R}");
4235  // fhDistanceHijingPhotonAddedSecondarySignal->SetXTitle("#it{p}_{T}");
4236  // outputContainer->Add(fhDistanceHijingPhotonAddedSecondarySignal) ;
4237  //
4238  // fhDistanceAddedPhotonHijingSecondary = new TH2F
4239  // ("hDistanceAddedPhotonHijingSecondary", "Distance added #gamma to secondary particle from hijing"
4240  // ,nptbins,ptmin,ptmax,100,0,0.4);
4241  // fhDistanceAddedPhotonHijingSecondary->SetYTitle("#it{R}");
4242  // fhDistanceAddedPhotonHijingSecondary->SetXTitle("#it{p}_{T}");
4243  // outputContainer->Add(fhDistanceAddedPhotonHijingSecondary) ;
4244 
4245 
4247  ("hDistance2AddedSignals", "Distance added signals"
4248  ,nptbins,ptmin,ptmax,100,0,0.4);
4249  fhDistance2AddedSignals->SetYTitle("#it{R}");
4250  fhDistance2AddedSignals->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4251  outputContainer->Add(fhDistance2AddedSignals) ;
4252 
4253  fhDistance2Hijing = new TH2F
4254  ("hDistance2Hijing", "Distance 2 hijing clusters"
4255  ,nptbins,ptmin,ptmax,100,0,0.4);
4256  fhDistance2Hijing->SetYTitle("#it{R}");
4257  fhDistance2Hijing->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4258  outputContainer->Add(fhDistance2Hijing) ;
4259 
4261  ("hDistanceAddedSignalsHijing", "Distance added signals to hijing"
4262  ,nptbins,ptmin,ptmax,100,0,0.4);
4263  fhDistanceAddedSignalsHijing->SetYTitle("#it{R}");
4264  fhDistanceAddedSignalsHijing->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4265  outputContainer->Add(fhDistanceAddedSignalsHijing) ;
4266  }
4267 
4269  {
4270 
4271  TString mcGenNames[] = {"","_MC_Pi0Merged","_MC_Pi0Decay","_MC_EtaDecay","_MC_PhotonOther","_MC_Electron","_MC_Other"};
4272  TString mcGenTitle[] = {"",",MC Pi0-Merged",",MC Pi0-Decay",", MC Eta-Decay",", MC Photon other sources",", MC Electron",", MC other sources"};
4273  for(Int_t igen = 0; igen < GetNCocktailGenNamesToCheck(); igen++)
4274  {
4275  TString add = "_MainGener_";
4276  if(igen==0) add = "";
4277  for(Int_t imc = 0; imc < fgkNGenTypes; imc++)
4278  {
4279  fhCleanGeneratorCluster[igen][imc] = new TH1F(Form("hCleanGeneratorCluster%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4280  Form("Number of selected clusters with contribution of %s generator%s, no overlap",
4281  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4282  nptbins,ptmin,ptmax);
4283  fhCleanGeneratorCluster[igen][imc]->SetYTitle("#it{counts}");
4284  fhCleanGeneratorCluster[igen][imc]->SetXTitle("#it{E} (GeV)");
4285  outputContainer->Add(fhCleanGeneratorCluster[igen][imc]) ;
4286 
4287  fhCleanGeneratorClusterEPrimRecoRatio[igen][imc] = new TH2F(Form("hCleanGeneratorClusterEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4288  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of %s generator%s, no overlap",
4289  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4290  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4291  fhCleanGeneratorClusterEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4292  fhCleanGeneratorClusterEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4293  outputContainer->Add(fhCleanGeneratorClusterEPrimRecoRatio[igen][imc]) ;
4294 
4295  fhCleanGeneratorClusterEPrimRecoDiff[igen][imc] = new TH2F(Form("hCleanGeneratorClusterEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4296  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of %s generator%s, no overlap",
4297  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4298  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4299  fhCleanGeneratorClusterEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4300  fhCleanGeneratorClusterEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4301  outputContainer->Add(fhCleanGeneratorClusterEPrimRecoDiff[igen][imc]) ;
4302 
4303  //
4304 
4305  fhMergeGeneratorCluster[igen][imc] = new TH1F(Form("hMergeGeneratorCluster%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4306  Form("Number of selected clusters with contribution of >=2 generators, main %s%s",
4307  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4308  nptbins,ptmin,ptmax);
4309  fhMergeGeneratorCluster[igen][imc]->SetYTitle("#it{counts}");
4310  fhMergeGeneratorCluster[igen][imc]->SetXTitle("#it{E} (GeV)");
4311  outputContainer->Add(fhMergeGeneratorCluster[igen][imc]) ;
4312 
4313  fhMergeGeneratorClusterEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4314  Form("#it{E}_{reco}/#it{E}_{gen}clusters with contribution of >=2 generators, main %s%s",
4315  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4316  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4317  fhMergeGeneratorClusterEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4318  fhMergeGeneratorClusterEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4319  outputContainer->Add(fhMergeGeneratorClusterEPrimRecoRatio[igen][imc]) ;
4320 
4321  fhMergeGeneratorClusterEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4322  Form("#it{E}_{reco}-#it{E}_{gen}clusters with contribution of >=2 generators, main %s%s",
4323  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4324  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4325  fhMergeGeneratorClusterEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4326  fhMergeGeneratorClusterEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4327  outputContainer->Add(fhMergeGeneratorClusterEPrimRecoDiff[igen][imc]) ;
4328 
4329  //if(GetCocktailGenNameToCheck(igen).Contains("ijing")) continue;
4330 
4331  //
4332 
4333  fhMergeGeneratorClusterNotHijingBkg[igen][imc] = new TH1F(Form("hMergeGeneratorClusterNotHijingBkg%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4334  Form("Number of selected clusters with contribution of >=2 generators, , none is HIJING, main %s%s",
4335  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4336  nptbins,ptmin,ptmax);
4337  fhMergeGeneratorClusterNotHijingBkg[igen][imc]->SetYTitle("#it{counts}");
4338  fhMergeGeneratorClusterNotHijingBkg[igen][imc]->SetXTitle("#it{E} (GeV)");
4339  outputContainer->Add(fhMergeGeneratorClusterNotHijingBkg[igen][imc]) ;
4340 
4341  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterNotHijingBkgEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4342  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of >=2 generators, , none is HIJING, main %s%s",
4343  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4344  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4345  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4346  fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4347  outputContainer->Add(fhMergeGeneratorClusterNotHijingBkgEPrimRecoRatio[igen][imc]) ;
4348 
4349  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterNotHijingBkgEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4350  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of >=2 generators, , none is HIJING, main %s%s",
4351  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4352  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4353  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4354  fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4355  outputContainer->Add(fhMergeGeneratorClusterNotHijingBkgEPrimRecoDiff[igen][imc]) ;
4356 
4357  //
4358 
4359  fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc] = new TH1F(Form("hMergeGeneratorClusterHijingAndOtherBkg%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4360  Form("Number of selected clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4361  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4362  nptbins,ptmin,ptmax);
4363  fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc]->SetYTitle("#it{counts}");
4364  fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc]->SetXTitle("#it{E} (GeV)");
4365  outputContainer->Add(fhMergeGeneratorClusterHijingAndOtherBkg[igen][imc]) ;
4366 
4367 
4368  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4369  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4370  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4371  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4372  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4373  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4374  outputContainer->Add(fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoRatio[igen][imc]) ;
4375 
4376 
4377  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4378  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4379  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4380  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4381  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4382  fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4383  outputContainer->Add(fhMergeGeneratorClusterHijingAndOtherBkgEPrimRecoDiff[igen][imc]) ;
4384 
4385  //
4386 
4387  fhMergeGeneratorClusterHijingBkg[igen][imc] = new TH1F(Form("hMergeGeneratorClusterHijingBkg%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4388  Form("Number of selected clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4389  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4390  nptbins,ptmin,ptmax);
4391  fhMergeGeneratorClusterHijingBkg[igen][imc]->SetYTitle("#it{counts}");
4392  fhMergeGeneratorClusterHijingBkg[igen][imc]->SetXTitle("#it{E} (GeV)");
4393  outputContainer->Add(fhMergeGeneratorClusterHijingBkg[igen][imc]) ;
4394 
4395  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingBkgEPrimRecoRatio%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4396  Form("#it{E}_{reco}/#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4397  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4398  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
4399  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
4400  fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4401  outputContainer->Add(fhMergeGeneratorClusterHijingBkgEPrimRecoRatio[igen][imc]) ;
4402 
4403  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc] = new TH2F(Form("hMergeGeneratorClusterHijingBkgEPrimRecoDiff%s%s%s",add.Data(),GetCocktailGenNameToCheck(igen).Data(), mcGenNames[imc].Data()),
4404  Form("#it{E}_{reco}-#it{E}_{gen} clusters with contribution of >=3 generators, none is HIJING, main %s%s",
4405  GetCocktailGenNameToCheck(igen).Data(), mcGenTitle[imc].Data()),
4406  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
4407  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
4408  fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc]->SetXTitle("#it{E}_{reco} (GeV)");
4409  outputContainer->Add(fhMergeGeneratorClusterHijingBkgEPrimRecoDiff[igen][imc]) ;
4410  }
4411  }
4412  }
4413 
4414  return outputContainer ;
4415 }
4416 
4417 //_______________________
4419 //_______________________
4421 {
4422  if ( GetCalorimeter() == kPHOS && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD() )
4423  AliFatal("!!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
4424  else if( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
4425  AliFatal("!!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
4426 
4427  // Trick to select primary photon particles in the analysis of pure MC input
4428  if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
4429 }
4430 
4431 //_________________________________
4433 //_________________________________
4435 {
4436  AddToHistogramsName("AnaPhoton_");
4437 
4438  fMinDist = 2.;
4439  fMinDist2 = 4.;
4440  fMinDist3 = 5.;
4441 
4442  fTimeCutMin =-1000000;
4443  fTimeCutMax = 1000000;
4444  fNCellsCut = 0;
4445 
4446  fRejectTrackMatch = kTRUE ;
4447 
4448  fNEBinCuts = 14;
4449  fEBinCuts[0] = 0.; fEBinCuts[1] = 0.3; fEBinCuts[2] = 0.5;
4450  fEBinCuts[3] = 1.; fEBinCuts[4] = 2. ; fEBinCuts[5] = 3. ;
4451  fEBinCuts[6] = 4.; fEBinCuts[7] = 5. ; fEBinCuts[8] = 7. ;
4452  fEBinCuts[9] = 9.; fEBinCuts[10]= 12.; fEBinCuts[11]= 15.;
4453  fEBinCuts[12]= 20.; fEBinCuts[13]= 50.; fEBinCuts[14]= 100.;
4454  for(Int_t i = fNEBinCuts; i < 15; i++) fEBinCuts[i] = 1000.;
4455 }
4456 
4457 //_______________________________________
4461 //_______________________________________
4463 {
4464  // Get the vertex
4465  Double_t v[3] = {0,0,0}; //vertex ;
4466  GetReader()->GetVertex(v);
4467 
4468  // Select the Calorimeter of the photon
4469  TObjArray * pl = 0x0;
4470  AliVCaloCells* cells = 0;
4471  if (GetCalorimeter() == kPHOS )
4472  {
4473  pl = GetPHOSClusters();
4474  cells = GetPHOSCells();
4475  }
4476  else if (GetCalorimeter() == kEMCAL)
4477  {
4478  pl = GetEMCALClusters();
4479  cells = GetEMCALCells();
4480  }
4481 
4482  if(!pl)
4483  {
4484  AliWarning(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
4485  return;
4486  }
4487 
4488  // Loop on raw clusters before filtering in the reader and fill control histogram
4489  if((GetReader()->GetEMCALClusterListName()=="" && GetCalorimeter()==kEMCAL) || GetCalorimeter()==kPHOS)
4490  {
4491  for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
4492  {
4493  AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
4494 
4495  if (GetCalorimeter() == kPHOS && clus->IsPHOS() && clus->E() > GetReader()->GetPHOSPtMin() )
4496  {
4497  clus->GetMomentum(fMomentum,GetVertex(0)) ;
4498 
4499  fhClusterCutsE [0]->Fill(fMomentum.E() , GetEventWeight());
4500  fhClusterCutsPt[0]->Fill(fMomentum.Pt(), GetEventWeight());
4501  }
4502  else if(GetCalorimeter() == kEMCAL && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin())
4503  {
4504  clus->GetMomentum(fMomentum,GetVertex(0)) ;
4505 
4506  fhClusterCutsE [0]->Fill(fMomentum.E(), GetEventWeight());
4507  fhClusterCutsPt[0]->Fill(fMomentum.Pt(), GetEventWeight());
4508  }
4509  }
4510  }
4511  else
4512  { // reclusterized
4513  TClonesArray * clusterList = 0;
4514 
4515  if(GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()))
4516  clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
4517  else if(GetReader()->GetOutputEvent())
4518  clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
4519 
4520  if(clusterList)
4521  {
4522  Int_t nclusters = clusterList->GetEntriesFast();
4523  for (Int_t iclus = 0; iclus < nclusters; iclus++)
4524  {
4525  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
4526 
4527  if(clus && clus->E() > GetReader()->GetEMCALPtMin())
4528  {
4529  clus->GetMomentum(fMomentum,GetVertex(0)) ;
4530 
4531  fhClusterCutsE [0]->Fill(clus->E() , GetEventWeight());
4532  fhClusterCutsPt[0]->Fill(fMomentum.Pt(), GetEventWeight());
4533  }
4534  }
4535  }
4536  }
4537 
4538  // Init arrays, variables, get number of clusters
4539  Int_t nCaloClusters = pl->GetEntriesFast();
4540 
4541  AliDebug(1,Form("Input %s cluster entries %d", GetCalorimeterString().Data(), nCaloClusters));
4542 
4543  //----------------------------------------------------
4544  // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
4545  //----------------------------------------------------
4546  // Loop on clusters
4547  for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
4548  {
4549  AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
4550  //printf("calo %d, %f\n",icalo,calo->E());
4551 
4552  //Get the index where the cluster comes, to retrieve the corresponding vertex
4553  Int_t evtIndex = 0 ;
4554  if (GetMixedEvent())
4555  {
4556  evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
4557  //Get the vertex and check it is not too large in z
4558  if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
4559  }
4560 
4561  //Cluster selection, not charged, with photon id and in fiducial cut
4563  {
4564  calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
4565  }//Assume that come from vertex in straight line
4566  else
4567  {
4568  Double_t vertex[]={0,0,0};
4569  calo->GetMomentum(fMomentum,vertex) ;
4570  }
4571 
4572  //-------------------------------------
4573  // Play with the MC stack if available
4574  //-------------------------------------
4575 
4576  // Check origin of the candidates
4577  Int_t tag = -1;
4578 
4579  if ( IsDataMC() )
4580  {
4581  tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(), GetMC(), pl); // check lost decays
4582 
4583  AliDebug(1,Form("Origin of candidate, bit map %d",tag));
4584  }// Work with stack also
4585 
4586  //-----------------------------
4587  // Cluster selection
4588  //-----------------------------
4589  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
4590  if(!ClusterSelected(calo,nMaxima,tag)) continue;
4591 
4592  //----------------------------
4593  // Create AOD for analysis
4594  //----------------------------
4595  AliAODPWG4Particle aodph = AliAODPWG4Particle(fMomentum);
4596 
4597  //...............................................
4598  // Set the indeces of the original caloclusters (MC, ID), and calorimeter
4599  Int_t label = calo->GetLabel();
4600  aodph.SetLabel(label);
4601  aodph.SetCaloLabel(calo->GetID(),-1);
4602  aodph.SetTag(tag);
4603  aodph.SetDetectorTag(GetCalorimeter());
4604  //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
4605 
4606  //...............................................
4607  // Set bad channel distance bit
4608  Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
4609  if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
4610  else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
4611  else aodph.SetDistToBad(0) ;
4612  //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
4613 
4614  //--------------------------------------------------------
4615  // Fill some shower shape histograms before PID is applied
4616  //--------------------------------------------------------
4617 
4618  Float_t maxCellFraction = 0;
4619  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo, maxCellFraction);
4620  if( absIdMax < 0 ) AliFatal("Wrong absID");
4621 
4622  Int_t largeTimeInCellCluster = kFALSE;
4623  FillShowerShapeHistograms(calo,tag,nMaxima,maxCellFraction,largeTimeInCellCluster);
4624  aodph.SetFiducialArea(largeTimeInCellCluster); // Temporary use of this container, FIXME
4625  //if(largeTimeInCellCluster > 1) printf("Set n cells large time %d, pt %2.2f\n",aodph.GetFiducialArea(),aodph.Pt());
4626 
4627  aodph.SetM02(calo->GetM02());
4628  aodph.SetM20(calo->GetM20());
4629  aodph.SetNLM(nMaxima);
4630 
4631  Float_t time = calo->GetTOF()*1e9;
4632  if(time > 400) time-=fConstantTimeShift; // in case of clusterizer running before (calibrate clusters not cells)
4633  aodph.SetTime(time);
4634 
4635  aodph.SetNCells(calo->GetNCells());
4636  Int_t nSM = GetModuleNumber(calo);
4637  aodph.SetSModNumber(nSM);
4638 
4639  Float_t en = fMomentum.E ();
4640  Float_t pt = fMomentum.Pt();
4641  Float_t eta = fMomentum.Eta();
4642  Float_t phi = GetPhi(fMomentum.Phi());
4643  Int_t ebin = -1;
4644  for(Int_t ie = 0; ie < fNEBinCuts; ie++)
4645  {
4646  if( en >= fEBinCuts[ie] && en < fEBinCuts[ie+1] ) ebin = ie;
4647  }
4648 
4649  Int_t icolAbs = -1, irowAbs = -1;
4651  {
4652  Float_t maxCellFraction = 0;
4653  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
4654 
4655  Int_t icol = -1, irow = -1, iRCU = -1;
4656  GetModuleNumberCellIndexesAbsCaloMap(absIdMax,GetCalorimeter(), icol, irow, iRCU, icolAbs, irowAbs);
4657 
4658  if(ebin>=0 && ebin < fNEBinCuts)
4659  {
4660  fhEBinClusterEtaPhi[ebin]->Fill(eta,phi,GetEventWeight()) ;
4661 
4662  fhEBinClusterColRow[ebin]->Fill(icolAbs,irowAbs,GetEventWeight()) ;
4663  }
4664  }
4665 
4666  //-------------------------------------
4667  // PID selection or bit setting
4668  //-------------------------------------
4669 
4670  //...............................................
4671  // Data, PID check on
4672  if(IsCaloPIDOn())
4673  {
4674  // Get most probable PID, 2 options check bayesian PID weights or redo PID
4675  // By default, redo PID
4676 
4677  aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));
4678 
4679  AliDebug(1,Form("PDG of identified particle %d",aodph.GetIdentifiedParticleType()));
4680 
4681  //If cluster does not pass pid, not photon, skip it.
4682  if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
4683  }
4684 
4685  //...............................................
4686  // Data, PID check off
4687  else
4688  {
4689  // Set PID bits for later selection (AliAnaPi0 for example)
4690  // GetIdentifiedParticleType already called in SetPIDBits.
4691 
4692  GetCaloPID()->SetPIDBits(calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
4693 
4694  AliDebug(1,"PID Bits set");
4695  }
4696 
4697  AliDebug(1,Form("Photon selection cuts passed: pT %3.2f, pdg %d",aodph.Pt(),aodph.GetIdentifiedParticleType()));
4698 
4699  fhClusterCutsE [9]->Fill(en, GetEventWeight());
4700  fhClusterCutsPt[9]->Fill(pt, GetEventWeight());
4701 
4702  //
4703  // Check local cluster activity around the current cluster
4704  //
4705  if(fStudyActivityNearCluster && en > 1.5) // 1.5 GeV cut used on Pb-Pb analysis
4706  ActivityNearCluster(icalo,en,eta,phi,tag,pl);
4707 
4708  //
4709  // Check if other generators contributed to the cluster
4710  //
4713 
4714 
4716  {
4717  if(ebin>=0 && ebin < fNEBinCuts)
4718  {
4719  fhEBinClusterEtaPhiPID[ebin]->Fill(eta,phi,GetEventWeight()) ;
4720 
4721  fhEBinClusterColRowPID[ebin]->Fill(icolAbs,irowAbs,GetEventWeight()) ;
4722  }
4723  }
4724 
4725  if(nSM < fNModules && nSM >=0)
4726  {
4727  fhEPhotonSM ->Fill(en, nSM, GetEventWeight());
4728  fhPtPhotonSM->Fill(pt, nSM, GetEventWeight());
4729  }
4730 
4731  // Few more control histograms for selected clusters
4732  fhNCellsE ->Fill(en, calo->GetNCells() , GetEventWeight());
4733  fhTimePt ->Fill(pt, time , GetEventWeight());
4734  fhNLocMax ->Fill(en, nMaxima , GetEventWeight());
4735 
4737  fhMaxCellDiffClusterE->Fill(en, maxCellFraction, GetEventWeight());
4738 
4739  if(cells)
4740  {
4741  for(