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