AliPhysics  f57202d (f57202d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaCalorimeterQA.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 is 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 <TObjArray.h>
18 #include <TParticle.h>
19 #include <TDatabasePDG.h>
20 #include <TH3F.h>
21 #include <TObjString.h>
22 
23 //---- AliRoot system ----
24 #include "AliAnaCalorimeterQA.h"
25 #include "AliCaloTrackReader.h"
26 #include "AliStack.h"
27 #include "AliVCaloCells.h"
28 #include "AliFiducialCut.h"
29 #include "AliVCluster.h"
30 #include "AliVTrack.h"
31 #include "AliVEvent.h"
32 #include "AliVEventHandler.h"
33 #include "AliAODMCParticle.h"
34 #include "AliMCAnalysisUtils.h"
35 
36 // --- Detectors ---
37 #include "AliPHOSGeoUtils.h"
38 #include "AliEMCALGeometry.h"
39 
41 ClassImp(AliAnaCalorimeterQA) ;
43 
44 //__________________________________________
46 //__________________________________________
47 AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
48 AliAnaCaloTrackCorrBaseClass(),
49 
50 // Switches
51 fFillAllCellTimeHisto(kTRUE),
52 fFillAllPosHisto(kFALSE), fFillAllPosHisto2(kFALSE),
53 fFillAllTH3(kFALSE),
54 fFillAllTMHisto(kTRUE),
55 fFillAllPi0Histo(kTRUE), fFillInvMassOpenAngle(kFALSE),
56 fFillPi0PairDiffTime(kFALSE), fFillInvMassInEMCALWithPHOSDCalAcc(kFALSE),
57 fFillEBinAcceptanceHisto(kFALSE),
58 fCorrelate(kTRUE), fStudyBadClusters(kFALSE),
59 fStudyClustersAsymmetry(kFALSE), fStudyExotic(kFALSE),
60 fStudyWeight(kFALSE), fStudyFECCorrelation(kFALSE),
61 fStudyM02Dependence (kFALSE),
62 
63 // Parameters and cuts
64 fNModules(12), fNRCU(2),
65 fNMaxCols(48), fNMaxRows(24),
66 fTimeCutMin(-10000), fTimeCutMax(10000),
67 fCellAmpMin(0),
68 fEMCALCellAmpMin(0), fPHOSCellAmpMin(0),
69 fEMCALClusterNCellMin(0), fPHOSClusterNCellMin(0),
70 fNEBinCuts(0),
71 
72 // Invariant mass
73 fInvMassMinECut(0), fInvMassMaxECut(0),
74 fInvMassMinM02Cut(0), fInvMassMaxM02Cut(0),
75 fInvMassMaxOpenAngle(0), fInvMassMaxTimeDifference(0),
76 
77 // Exotic
78 fExoNECrossCuts(0), fExoECrossCuts(),
79 fExoNDTimeCuts(0), fExoDTimeCuts(),
80 
81 fClusterMomentum(), fClusterMomentum2(),
82 fPrimaryMomentum(),
83 fConstantTimeShift(0),
84 // Histograms
85 fhE(0), fhPt(0),
86 fhPhi(0), fhEta(0),
87 fhEtaPhi(0), fhEtaPhiE(0),
88 fhECharged(0), fhPtCharged(0),
89 fhPhiCharged(0), fhEtaCharged(0),
90 fhEtaPhiCharged(0), fhEtaPhiECharged(0),
91 
92 // Invariant mass
93 fhIM(0), fhIMSame(0), fhIMDiff(0),
94 fhIMDCAL(0), fhIMDCALSame(0), fhIMDCALDiff(0),
95 fhIMDCALPHOS(0), fhIMDCALPHOSSame(0),
96 fhIMEMCALPHOS(0), fhIMEMCALPHOSSame(0),
97 fhAsym(0),
98 fhOpAngle(0), fhIMvsOpAngle(0),
99 fhNCellsPerCluster(0), fhNCellsPerClusterNoCut(0),
100 fhNCellsPerClusterWithWeight(0), fhNCellsPerClusterRatioWithWeight(0),
101 fhNClusters(0),
102 
103 // Timing
104 fhClusterTimeEnergy(0), fhCellTimeSpreadRespectToCellMax(0),
105 fhCellIdCellLargeTimeSpread(0), fhClusterPairDiffTimeE(0), fhClusterPairDiffTimeESameMod(0),
106 fhClusterMaxCellCloseCellRatio(0), fhClusterMaxCellCloseCellDiff(0),
107 fhClusterMaxCellDiff(0), fhClusterMaxCellDiffNoCut(0),
108 //fhClusterMaxCellDiffAverageTime(0), fhClusterMaxCellDiffWeightedTime(0),
109 fhClusterMaxCellECross(0),
110 fhLambda0(0), fhLambda1(0),
111 fhLambda0MaxFECPair(0), fhLambda1MaxFECPair(0),
112 fhLambda0MaxFECOdd(0), fhLambda1MaxFECOdd(0),
113 fhLambda0Max2FECPair(0), fhLambda1Max2FECPair(0),
114 fhLambda0Max2FECOdd(0), fhLambda1Max2FECOdd(0),
115 fhLambda0MaxFECPairLargeNCells(0), fhLambda1MaxFECPairLargeNCells(0),
116 fhLambda0MaxFECOddLargeNCells(0), fhLambda1MaxFECOddLargeNCells(0),
117 fhLambda0Max2FECPairLargeNCells(0), fhLambda1Max2FECPairLargeNCells(0),
118 fhLambda0Max2FECOddLargeNCells(0), fhLambda1Max2FECOddLargeNCells(0),
119 
120 // fhDispersion(0),
121 fhEtaPhiFECCorrControl(0),
122 // bad clusters
123 fhBadClusterEnergy(0), fhBadClusterTimeEnergy(0), fhBadClusterEtaPhi(0),
124 fhBadClusterPairDiffTimeE(0), fhBadCellTimeSpreadRespectToCellMax(0),
125 fhBadClusterMaxCellCloseCellRatio(0), fhBadClusterMaxCellCloseCellDiff(0), fhBadClusterMaxCellDiff(0),
126 fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
127 fhBadClusterMaxCellECross(0),
128 fhBadClusterLambda0(0), fhBadClusterLambda1(0),
129 fhBadClusterDeltaIEtaDeltaIPhiE0(0), fhBadClusterDeltaIEtaDeltaIPhiE2(0),
130 fhBadClusterDeltaIEtaDeltaIPhiE6(0), fhBadClusterDeltaIA(0),
131 
132 // Position
133 fhRNCells(0), fhXNCells(0),
134 fhYNCells(0), fhZNCells(0),
135 fhRE(0), fhXE(0),
136 fhYE(0), fhZE(0),
137 fhXYZ(0),
138 fhRCellE(0), fhXCellE(0),
139 fhYCellE(0), fhZCellE(0),
140 fhXYZCell(0),
141 fhDeltaCellClusterRNCells(0), fhDeltaCellClusterXNCells(0),
142 fhDeltaCellClusterYNCells(0), fhDeltaCellClusterZNCells(0),
143 fhDeltaCellClusterRE(0), fhDeltaCellClusterXE(0),
144 fhDeltaCellClusterYE(0), fhDeltaCellClusterZE(0),
145 
146 // Cells
147 fhNCells(0), fhNCellsCutAmpMin(0),
148 fhAmplitude(0), fhAmpId(0),
149 fhEtaPhiAmpCell(0), fhEtaPhiCell(0),
150 fhTime(0), //fhTimeVz(0),
151 fhTimeId(0), fhTimeAmp(0),
152 fhAmpIdLowGain(0), fhTimeIdLowGain(0), fhTimeAmpLowGain(0),
153 
154 fhCellECross(0),
155 
156 fhEMCALPHOSCorrNClusters(0), fhEMCALPHOSCorrEClusters(0),
157 fhEMCALPHOSCorrNCells(0), fhEMCALPHOSCorrECells(0),
158 fhEMCALDCALCorrNClusters(0), fhEMCALDCALCorrEClusters(0),
159 fhEMCALDCALCorrNCells(0), fhEMCALDCALCorrECells(0),
160 fhDCALPHOSCorrNClusters(0), fhDCALPHOSCorrEClusters(0),
161 fhDCALPHOSCorrNCells(0), fhDCALPHOSCorrECells(0),
162 fhCaloV0SCorrNClusters(0), fhCaloV0SCorrEClusters(0),
163 fhCaloV0SCorrNCells(0), fhCaloV0SCorrECells(0),
164 fhCaloV0MCorrNClusters(0), fhCaloV0MCorrEClusters(0),
165 fhCaloV0MCorrNCells(0), fhCaloV0MCorrECells(0),
166 fhCaloTrackMCorrNClusters(0), fhCaloTrackMCorrEClusters(0),
167 fhCaloTrackMCorrNCells(0), fhCaloTrackMCorrECells(0),
168 fhCaloCenNClusters(0), fhCaloCenEClusters(0),
169 fhCaloCenNCells(0), fhCaloCenECells(0),
170 fhCaloEvPNClusters(0), fhCaloEvPEClusters(0),
171 fhCaloEvPNCells(0), fhCaloEvPECells(0),
172 
173 // Super-Module dependent histograms
174 fhEMod(0), fhAmpMod(0),
175 fhEWeirdMod(0), fhAmpWeirdMod(0),
176 fhTimeMod(0),
177 fhNClustersMod(0), fhNCellsMod(0),
178 fhNCellsSumAmpPerMod(0), fhNClustersSumEnergyPerMod(0),
179 fhNCellsPerClusterMod(0), fhNCellsPerClusterModNoCut(0),
180 fhNCellsPerClusterWeirdMod(0), fhNCellsPerClusterWeirdModNoCut(0),
181 
182 fhGridCells(0), fhGridCellsE(0), fhGridCellsTime(0),
183 fhGridCellsLowGain(0), fhGridCellsELowGain(0), fhGridCellsTimeLowGain(0),
184 fhTimeAmpPerRCU(0), fhIMMod(0),
185 
186 // Weight studies
187 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
188 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
189 fhECellTotalRatio(0), fhECellTotalLogRatio(0),
190 fhECellTotalRatioMod(0), fhECellTotalLogRatioMod(0),
191 
192 fhExoL0ECross(0), fhExoL1ECross(0),
193 
194 // MC and reco
195 fhRecoMCE(), fhRecoMCPhi(), fhRecoMCEta(),
196 fhRecoMCDeltaE(), fhRecoMCRatioE(),
197 fhRecoMCDeltaPhi(), fhRecoMCDeltaEta(),
198 
199 // MC only
200 fhGenMCE(), fhGenMCPt(), fhGenMCEtaPhi(),
201 fhGenMCAccE(), fhGenMCAccPt(), fhGenMCAccEtaPhi(),
202 
203 // Matched MC
204 fhEMVxyz(0), fhEMR(0),
205 fhHaVxyz(0), fhHaR(0),
206 fh1EOverP(0), fh2dR(0),
207 fh2EledEdx(0), fh2MatchdEdx(0),
208 fh1EOverPR02(0), fh1EleEOverP(0),
209 fhMCEle1EOverP(0), fhMCEle1dR(0), fhMCEle2MatchdEdx(0),
210 fhMCChHad1EOverP(0), fhMCChHad1dR(0), fhMCChHad2MatchdEdx(0),
211 fhMCNeutral1EOverP(0), fhMCNeutral1dR(0), fhMCNeutral2MatchdEdx(0),
212 fhMCEle1EOverPR02(0), fhMCChHad1EOverPR02(0), fhMCNeutral1EOverPR02(0),
213 fhMCEle1EleEOverP(0), fhMCChHad1EleEOverP(0), fhMCNeutral1EleEOverP(0),
214 fhTrackMatchedDEtaNeg(0), fhTrackMatchedDPhiNeg(0), fhTrackMatchedDEtaDPhiNeg(0),
215 fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), fhTrackMatchedDEtaDPhiPos(0),
216 fhTrackMatchedDEtaNegMod(0), fhTrackMatchedDPhiNegMod(0),
217 fhTrackMatchedDEtaPosMod(0), fhTrackMatchedDPhiPosMod(0),
218 fhClusterTimeEnergyM02(0), fhCellTimeSpreadRespectToCellMaxM02(0),
219 fhClusterMaxCellCloseCellRatioM02(0), fhClusterMaxCellCloseCellDiffM02(0),
220 fhClusterMaxCellDiffM02(0), fhClusterMaxCellECrossM02(0), fhNCellsPerClusterM02(0)
221 {
222  // Weight studies
223  for(Int_t i =0; i < 12; i++)
224  {
225  for(Int_t j = 0; j < 4; j++)
226  {
227  for(Int_t k = 0; k < 3; k++)
228  {
229  fhLambda0ForW0AndCellCuts [i][j][k] = 0;
230 // fhLambda1ForW0AndCellCuts [i][j][k] = 0;
231  fhLambda0ForW0AndCellCutsEta0[i][j][k] = 0;
232  }
233  }
234  for(Int_t j = 0; j < 5; j++)
235  {
236  fhLambda0ForW0MC[i][j] = 0;
237 // fhLambda1ForW0MC[i][j] = 0;
238  }
239  }
240 
241  // Cluster size
242  fhDeltaIEtaDeltaIPhiE0[0] = 0 ; fhDeltaIEtaDeltaIPhiE2[0] = 0; fhDeltaIEtaDeltaIPhiE6[0] = 0;
243  fhDeltaIEtaDeltaIPhiE0[1] = 0 ; fhDeltaIEtaDeltaIPhiE2[1] = 0; fhDeltaIEtaDeltaIPhiE6[1] = 0;
244  fhDeltaIA[0] = 0 ; fhDeltaIAL0[0] = 0; fhDeltaIAL1[0] = 0;
245  fhDeltaIA[1] = 0 ; fhDeltaIAL0[1] = 0; fhDeltaIAL1[1] = 0;
246  fhDeltaIANCells[0] = 0 ; fhDeltaIANCells[1] = 0;
247  fhDeltaIAMC[0] = 0 ; fhDeltaIAMC[1] = 0;
248  fhDeltaIAMC[2] = 0 ; fhDeltaIAMC[3] = 0;
249 
250  // Exotic
251  for (Int_t ie = 0; ie < 10 ; ie++)
252  {
253  fhExoDTime[ie] = 0;
254  for (Int_t idt = 0; idt < 5 ; idt++)
255  {
256  fhExoNCell [ie][idt] = 0;
257  fhExoL0 [ie][idt] = 0;
258  fhExoL1 [ie][idt] = 0;
259  fhExoECross [ie][idt] = 0;
260  fhExoTime [ie][idt] = 0;
261  fhExoL0NCell [ie][idt] = 0;
262  fhExoL1NCell [ie][idt] = 0;
263  }
264  }
265 
266  // MC
267 
268  for(Int_t i = 0; i < 7; i++)
269  {
270  fhRecoMCE[i][0] = 0; fhRecoMCE[i][1] = 0;
271  fhRecoMCPhi[i][0] = 0; fhRecoMCPhi[i][1] = 0;
272  fhRecoMCEta[i][0] = 0; fhRecoMCEta[i][1] = 0;
273  fhRecoMCDeltaE[i][0] = 0; fhRecoMCDeltaE[i][1] = 0;
274  fhRecoMCRatioE[i][0] = 0; fhRecoMCRatioE[i][1] = 0;
275  fhRecoMCDeltaPhi[i][0] = 0; fhRecoMCDeltaPhi[i][1] = 0;
276  fhRecoMCDeltaEta[i][0] = 0; fhRecoMCDeltaEta[i][1] = 0;
277  }
278 
279  for(Int_t i = 0; i < 4; i++)
280  {
281  fhGenMCE[i] = 0;
282  fhGenMCPt[i] = 0;
283  fhGenMCEtaPhi[i] = 0;
284  fhGenMCAccE[i] = 0;
285  fhGenMCAccPt[i] = 0;
286  fhGenMCAccEtaPhi[i] = 0;
287  }
288 
289  for(Int_t i = 0; i < 14; i++)
290  {
291  fhEBinClusterEtaPhi[i] = 0 ;
292  fhEBinClusterColRow[i] = 0 ;
293  fhEBinCellColRow [i] = 0 ;
294  }
295 
296  for(Int_t bc = 0; bc < 4; bc++) fhTimePerSMPerBC[bc] = 0 ;
297 
298  for(Int_t i = 0; i < 6; i++)
299  {
300  fhLambda0FECCorrel [i] = 0 ;
301  fhLambda1FECCorrel [i] = 0 ;
302  fhLambda0FECCorrelLargeNCells [i] = 0 ;
303  fhLambda1FECCorrelLargeNCells [i] = 0 ;
304  }
305  for(Int_t i = 0; i < 4; i++)
306  {
307  fhLambda0MaxFECCorrel[i] = 0 ;
308  fhLambda1MaxFECCorrel[i] = 0 ;
309  fhLambda0MaxFECCorrelLargeNCells[i] = 0 ;
310  fhLambda1MaxFECCorrelLargeNCells[i] = 0 ;
311  }
312 
313  InitParameters();
314 }
315 
316 //______________________________________________________________________________________________________________________
326 //______________________________________________________________________________________________________________________
327 void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
328  Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
329  Double_t tmax)
330 {
331  // printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
332  // GetReader()->GetEventNumber(), GetCalorimeterString().Data(),
333  // clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
334 
335  Double_t tof = clus->GetTOF()*1.e9;
336  if(tof > 400) tof-=fConstantTimeShift;
337 
338  Float_t energy = clus->E();
339 
340  fhBadClusterEnergy ->Fill(energy, GetEventWeight());
341  fhBadClusterTimeEnergy ->Fill(energy, tof , GetEventWeight());
342  fhBadClusterMaxCellDiff ->Fill(energy, maxCellFraction, GetEventWeight());
343  fhBadClusterMaxCellECross->Fill(energy, eCrossFrac , GetEventWeight());
344 
345  Float_t phi = fClusterMomentum.Phi();
346  if(phi < 0) phi += TMath::TwoPi();
347 
348  if(energy > 0.5) fhBadClusterEtaPhi->Fill(fClusterMomentum.Eta(), phi, GetEventWeight());
349 
350  fhBadClusterLambda0->Fill(energy, clus->GetM02());
351  fhBadClusterLambda1->Fill(energy, clus->GetM20());
352 
353  if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
354 
355  // Clusters in event time difference bad minus good
356 
357  for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
358  {
359  AliVCluster* clus2 = (AliVCluster*) caloClusters->At(iclus2);
360 
361  if(clus->GetID() == clus2->GetID()) continue;
362 
363  Float_t maxCellFraction2 = 0.;
364  Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
365 
366  if(IsGoodCluster(absIdMax2, clus->GetM02(), clus->GetNCells(), cells) && clus2->GetM02() > 0.1 )
367  {
368  Double_t tof2 = clus2->GetTOF()*1.e9;
369  if(tof2>400) tof2-=fConstantTimeShift;
370 
371  fhBadClusterPairDiffTimeE ->Fill(clus->E(), (tof-tof2), GetEventWeight());
372  }
373  } // loop
374 
375  // Max cell compared to other cells in cluster
376  if(fFillAllCellTimeHisto)
377  {
378  // Get some time averages
379  Double_t timeAverages[2] = {0.,0.};
380  CalculateAverageTime(clus, cells, timeAverages);
381 
382  fhBadClusterMaxCellDiffAverageTime ->Fill(clus->E(), tmax-timeAverages[0], GetEventWeight());
383  fhBadClusterMaxCellDiffWeightedTime->Fill(clus->E(), tmax-timeAverages[1], GetEventWeight());
384  }
385 
386  for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
387  {
388  Int_t absId = clus->GetCellsAbsId()[ipos];
389  if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
390  {
391  Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
392 
393  fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(), frac, GetEventWeight());
394  fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(), cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId), GetEventWeight());
395 
396  if(fFillAllCellTimeHisto)
397  {
398  Double_t time = cells->GetCellTime(absId);
399  GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
400 
401  Float_t diff = (tmax-(time*1e9-fConstantTimeShift));
402  fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff, GetEventWeight());
403 
404  }
405  } // Not max
406  } // loop
407 }
408 
409 //______________________________________________________________________
414 //______________________________________________________________________
415 void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus,
416  AliVCaloCells* cells,
417  Double_t timeAverages[2])
418 {
419  // First recalculate energy in case non linearity was applied
420  Float_t energy = 0;
421  Float_t ampMax = 0, amp = 0;
422 //Int_t absIdMax = -1;
423 
424  for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
425  {
426  Int_t id = clus->GetCellsAbsId()[ipos];
427 
428  // Recalibrate cell energy if needed
429  amp = cells->GetCellAmplitude(id);
430  GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
431 
432  energy += amp;
433 
434  if(amp> ampMax)
435  {
436  ampMax = amp;
437 // absIdMax = id;
438  }
439  } // energy loop
440 
441  // Calculate average time of cells in cluster and weighted average
442  Double_t aTime = 0;
443  Double_t wTime = 0;
444  Float_t wTot = 0;
445  Double_t time = 0;
446  Int_t id =-1;
447  Double_t w = 0;
448  Int_t ncells = clus->GetNCells();
449 
450  for (Int_t ipos = 0; ipos < ncells; ipos++)
451  {
452  id = clus ->GetCellsAbsId()[ipos];
453  amp = cells->GetCellAmplitude(id);
454  time = cells->GetCellTime(id);
455 
456  // Recalibrate energy and time
457  GetCaloUtils()->RecalibrateCellAmplitude(amp , GetCalorimeter(), id);
458  GetCaloUtils()->RecalibrateCellTime (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
459 
460  w = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
461  aTime += time*1e9;
462  wTime += time*1e9 * w;
463  wTot += w;
464  }
465 
466  if(ncells > 0) aTime /= ncells;
467  else aTime = 0;
468 
469  if(wTot > 0) wTime /= wTot;
470  else wTime = 0;
471 
472  timeAverages[0] = aTime;
473  timeAverages[1] = wTime;
474 }
475 
476 //____________________________________________________________
478 //____________________________________________________________
479 void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
480 {
481  Int_t ncells = cells->GetNumberOfCells();
482  if( ncells > 0 ) fhNCells->Fill(ncells, GetEventWeight()) ; // Not ok for PHOS with CPV
483 
484  Int_t ncellsCut = 0;
485  Float_t ecellsCut = 0;
486 
487  AliDebug(1,Form("%s cell entries %d", GetCalorimeterString().Data(), ncells));
488 
489  // Init arrays and used variables
490  Int_t *nCellsInModule = new Int_t [fNModules];
491  Float_t *eCellsInModule = new Float_t[fNModules];
492 
493  for(Int_t imod = 0; imod < fNModules; imod++ )
494  {
495  nCellsInModule[imod] = 0 ;
496  eCellsInModule[imod] = 0.;
497  }
498 
499  Int_t icol = -1, icolAbs = -1;
500  Int_t irow = -1, irowAbs = -1;
501  Int_t iRCU = -1;
502  Float_t amp = 0.;
503  Double_t time = 0.;
504  Int_t id = -1;
505  Bool_t highG = kFALSE;
506  Float_t recalF = 1.;
507  Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
508 
509  for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
510  {
511  if ( cells->GetCellNumber(iCell) < 0 ) continue; // CPV
512 
513  AliDebug(2,Form("Cell : amp %f, absId %d", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell)));
514 
515  Int_t nModule = GetModuleNumberCellIndexesAbsCaloMap(cells->GetCellNumber(iCell),GetCalorimeter(),
516  icol , irow, iRCU,
517  icolAbs, irowAbs );
518 
519  AliDebug(2,Form("\t module %d, column %d (%d), row %d (%d)", nModule,icolAbs,icol,irowAbs,irow));
520 
521  if(nModule < fNModules)
522  {
523  //Check if the cell is a bad channel
524  if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
525  {
526  if(GetCalorimeter()==kEMCAL)
527  {
528  if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
529  }
530  else
531  {
532  if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
533  }
534  } // use bad channel map
535 
536  amp = cells->GetAmplitude(iCell)*recalF;
537  time = cells->GetTime(iCell);
538  id = cells->GetCellNumber(iCell);
539  highG = cells->GetCellHighGain(id);
540  if(IsDataMC()) highG = kTRUE; // MC does not distinguish High and Low, put them all in high
541 
542  // Amplitude recalibration if set
543  GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
544 
545  // Time recalibration if set
546  GetCaloUtils()->RecalibrateCellTime (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
547 
548  // Transform time to ns
549  time *= 1.0e9;
550  time-=fConstantTimeShift;
551 
552  if(time < fTimeCutMin || time > fTimeCutMax)
553  {
554  AliDebug(1,Form("Remove cell with Time %f",time));
555  continue;
556  }
557 
558  // Remove exotic cells, defined only for EMCAL
559  if(GetCalorimeter()==kEMCAL &&
560  GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
561 
562  fhAmplitude->Fill(amp, GetEventWeight());
563  fhAmpId ->Fill(amp, id , GetEventWeight());
564  fhAmpMod ->Fill(amp, nModule, GetEventWeight());
565  fhAmpWeirdMod->Fill(amp, nModule, GetEventWeight());
566  if(!highG) fhAmpIdLowGain->Fill(amp, id, GetEventWeight());
567 
568  if(fFillEBinAcceptanceHisto)
569  {
570  for(Int_t ie = 0; ie < fNEBinCuts; ie++)
571  {
572  if( amp >= fEBinCuts[ie] && amp < fEBinCuts[ie+1] )
573  {
574  fhEBinCellColRow[ie]->Fill(icolAbs,irowAbs,GetEventWeight()) ;
575  }
576  }
577  }
578 
579  // E cross for exotic cells
580  if(amp > 0.05)
581  {
582  fhCellECross->Fill(amp, 1-GetECross(id,cells)/amp, GetEventWeight());
583  ecellsCut+=amp ;
584  if(fStudyWeight) eCellsInModule[nModule]+=amp ;
585  }
586 
587  if ( amp > fCellAmpMin )
588  {
589  ncellsCut++ ;
590  nCellsInModule[nModule]++ ;
591 
592  if(!fStudyWeight) eCellsInModule[nModule]+=amp ;
593 
594  fhGridCells ->Fill(icolAbs, irowAbs, GetEventWeight());
595  fhGridCellsE->Fill(icolAbs, irowAbs, amp );
596 
597  if(!highG)
598  {
599  fhGridCellsLowGain ->Fill(icolAbs, irowAbs, GetEventWeight());
600  fhGridCellsELowGain->Fill(icolAbs, irowAbs, amp );
601  }
602 
603  if(fFillAllCellTimeHisto)
604  {
605  //printf("%s: time %g\n",GetCalorimeterString().Data(), time);
606 
607 // Double_t v[3] = {0,0,0}; //vertex ;
608 // GetReader()->GetVertex(v);
609 // if(amp > 0.5) fhTimeVz ->Fill(TMath::Abs(v[2]), time, GetEventWeight());
610 
611  fhTime ->Fill(time, GetEventWeight());
612  fhTimeId ->Fill(time, id , GetEventWeight());
613  fhTimeAmp->Fill(amp , time, GetEventWeight());
614 
615  Int_t bc = (GetReader()->GetInputEvent()->GetBunchCrossNumber())%4;
616  fhTimePerSMPerBC[bc]->Fill(time, nModule, GetEventWeight());
617 
618  fhGridCellsTime->Fill(icolAbs, irowAbs, time);
619  if(!highG) fhGridCellsTimeLowGain->Fill(icolAbs, irowAbs, time);
620 
621  fhTimeMod->Fill(time, nModule, GetEventWeight());
622  fhTimeAmpPerRCU[nModule*fNRCU+iRCU]->Fill(amp, time, GetEventWeight());
623 
624  if(!highG)
625  {
626  fhTimeIdLowGain ->Fill(time, id , GetEventWeight());
627  fhTimeAmpLowGain->Fill(amp , time, GetEventWeight());
628  }
629  }
630  }
631 
632  // Get Eta-Phi position of Cell
633  if(fFillAllPosHisto)
634  {
635  if ( GetCalorimeter() == kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet() )
636  {
637  Float_t celleta = 0.;
638  Float_t cellphi = 0.;
639  GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
640 
641  if ( cellphi < 0 ) cellphi+=TMath::TwoPi();
642 
643  if(fFillAllTH3)
644  fhEtaPhiAmpCell->Fill(celleta, cellphi, amp, GetEventWeight());
645  else
646  fhEtaPhiCell ->Fill(celleta, cellphi, GetEventWeight());
647 
648  Double_t cellpos[] = {0, 0, 0};
649  GetEMCALGeometry()->GetGlobal(id, cellpos);
650 
651  fhXCellE->Fill(cellpos[0], amp, GetEventWeight()) ;
652  fhYCellE->Fill(cellpos[1], amp, GetEventWeight()) ;
653  fhZCellE->Fill(cellpos[2], amp, GetEventWeight()) ;
654 
655  Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
656  fhRCellE ->Fill(rcell, amp, GetEventWeight()) ;
657 
658  fhXYZCell->Fill(cellpos[0], cellpos[1], cellpos[2], GetEventWeight()) ;
659  } // EMCAL Cells
660  else if ( GetCalorimeter() == kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet() )
661  {
662  TVector3 xyz;
663  Int_t relId[4], module;
664  Float_t xCell, zCell;
665 
666  GetPHOSGeometry()->AbsToRelNumbering(id,relId);
667  module = relId[0];
668  GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
669  GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
670 
671  Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
672 
673  fhXCellE ->Fill(xyz.X(), amp, GetEventWeight()) ;
674  fhYCellE ->Fill(xyz.Y(), amp, GetEventWeight()) ;
675  fhZCellE ->Fill(xyz.Z(), amp, GetEventWeight()) ;
676  fhRCellE ->Fill(rcell , amp, GetEventWeight()) ;
677 
678  fhXYZCell->Fill(xyz.X(), xyz.Y(), xyz.Z(), GetEventWeight()) ;
679  } // PHOS cells
680  } // Fill cell position histograms
681 
682  } // N modules
683  } // Cell loop
684 
685  // Fill the cells after the cut on min amplitude and bad/exotic channels
686  if( ncellsCut > 0 ) fhNCellsCutAmpMin->Fill(ncellsCut, GetEventWeight()) ;
687 
688  // Number of cells per module
689  for(Int_t imod = 0; imod < fNModules; imod++ )
690  {
691  AliDebug(1,Form("Module %d, calo %s, N cells %d, sum Amp %f", imod, GetCalorimeterString().Data(), nCellsInModule[imod], eCellsInModule[imod]));
692 
693  fhNCellsMod ->Fill(nCellsInModule[imod], imod, GetEventWeight()) ;
694  fhSumCellsAmpMod->Fill(eCellsInModule[imod], imod, GetEventWeight()) ;
695 
696  fhNCellsSumAmpPerMod[imod]->Fill(eCellsInModule[imod], nCellsInModule[imod], GetEventWeight());
697  }
698 
699  // Check energy distribution in calorimeter for selected cells
700  if(fStudyWeight)
701  {
702  for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
703  {
704  if ( cells->GetCellNumber(iCell) < 0 ) continue; // CPV
705 
706  AliDebug(2,Form("Cell : amp %f, absId %d", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell)));
707 
708  Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
709 
710  AliDebug(2,Form("\t module %d, column %d, row %d", nModule,icol,irow));
711 
712  if(nModule < fNModules)
713  {
714  //Check if the cell is a bad channel
715  if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
716  {
717  if(GetCalorimeter()==kEMCAL)
718  {
719  if(GetCaloUtils()->GetEMCALChannelStatus(nModule, icol, irow)) continue;
720  }
721  else
722  {
723  if(GetCaloUtils()->GetPHOSChannelStatus (nModule, icol, irow) ) continue;
724  }
725  } // use bad channel map
726 
727  amp = cells->GetAmplitude(iCell)*recalF;
728  time = cells->GetTime(iCell);
729  id = cells->GetCellNumber(iCell);
730 
731  // Amplitude recalibration if set
732  GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
733 
734  // Time recalibration if set
735  GetCaloUtils()->RecalibrateCellTime (time, GetCalorimeter(), id,
736  GetReader()->GetInputEvent()->GetBunchCrossNumber());
737 
738  // Transform time to ns
739  time *= 1.0e9;
740  time -= fConstantTimeShift;
741 
742  if(time < fTimeCutMin || time > fTimeCutMax)
743  {
744  AliDebug(1,Form("Remove cell with Time %f",time));
745  continue;
746  }
747 
748  // Remove exotic cells, defined only for EMCAL
749  if(GetCalorimeter()==kEMCAL &&
750  GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
751 
752  // E cross for exotic cells
753  if(amp > 0.05)
754  {
755  if(ecellsCut > 0)
756  {
757  Float_t ratio = amp/ecellsCut;
758  fhECellTotalRatio ->Fill(ecellsCut, ratio , GetEventWeight());
759  fhECellTotalLogRatio ->Fill(ecellsCut,TMath::Log(ratio), GetEventWeight());
760  }
761 
762  if(eCellsInModule[nModule] > 0)
763  {
764  Float_t ratioMod = amp/eCellsInModule[nModule];
765  fhECellTotalRatioMod [nModule]->Fill(eCellsInModule[nModule], ratioMod , GetEventWeight());
766  fhECellTotalLogRatioMod[nModule]->Fill(eCellsInModule[nModule],TMath::Log(ratioMod), GetEventWeight());
767  }
768  } // amp > 0.5
769  } // nMod > 0 < Max
770  } // cell loop
771  } // weight studies
772 
773  delete [] nCellsInModule;
774  delete [] eCellsInModule;
775 }
776 
777 //__________________________________________________________________________
779 //__________________________________________________________________________
780 void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
781 {
782  Int_t nCaloCellsPerCluster = clus->GetNCells();
783 
784  UShort_t * indexList = clus->GetCellsAbsId();
785 
786  Float_t pos[3];
787  clus->GetPosition(pos);
788 
789  Float_t clEnergy = clus->E();
790 
791  // Loop on cluster cells
792  for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
793  {
794  // printf("Index %d\n",ipos);
795  Int_t absId = indexList[ipos];
796 
797  //Get position of cell compare to cluster
798 
799  if ( GetCalorimeter() == kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet() )
800  {
801  Double_t cellpos[] = {0, 0, 0};
802  GetEMCALGeometry()->GetGlobal(absId, cellpos);
803 
804  fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0], nCaloCellsPerCluster, GetEventWeight()) ;
805  fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1], nCaloCellsPerCluster, GetEventWeight()) ;
806  fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2], nCaloCellsPerCluster, GetEventWeight()) ;
807 
808  fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0], clEnergy, GetEventWeight()) ;
809  fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1], clEnergy, GetEventWeight()) ;
810  fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2], clEnergy, GetEventWeight()) ;
811 
812  Float_t r = TMath::Sqrt(pos[0] *pos[0] + pos[1] * pos[1] );
813  Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
814 
815  fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster, GetEventWeight()) ;
816  fhDeltaCellClusterRE ->Fill(r-rcell, clEnergy , GetEventWeight()) ;
817  } // EMCAL and its matrices are available
818  else if ( GetCalorimeter() == kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet() )
819  {
820  TVector3 xyz;
821  Int_t relId[4], module;
822  Float_t xCell, zCell;
823 
824  GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
825  module = relId[0];
826  GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
827  GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
828 
829  fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(), nCaloCellsPerCluster, GetEventWeight()) ;
830  fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(), nCaloCellsPerCluster, GetEventWeight()) ;
831  fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(), nCaloCellsPerCluster, GetEventWeight()) ;
832 
833  fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(), clEnergy, GetEventWeight()) ;
834  fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(), clEnergy, GetEventWeight()) ;
835  fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(), clEnergy, GetEventWeight()) ;
836 
837  Float_t r = TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1] );
838  Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
839 
840  fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster, GetEventWeight()) ;
841  fhDeltaCellClusterRE ->Fill(r-rcell, clEnergy , GetEventWeight()) ;
842  } // PHOS and its matrices are available
843  } // cluster cell loop
844 }
845 
846 
847 //___________________________________________________
855 //___________________________________________________
856 void AliAnaCalorimeterQA::ChannelCorrelationInFEC(AliVCluster* clus, AliVCaloCells* cells,
857  Bool_t matched, Int_t absIdMax) const
858 {
859  // Clean the sample
860 
861  // select neutral
862  if ( matched ) return;
863 
864  // away from dead region
865  if ( clus->GetDistanceToBadChannel() < 5 ) return ;
866 
867  // in center of SM
868  Int_t etaRegion = -1, phiRegion = -1;
869  GetCaloUtils()->GetEMCALSubregion(clus,GetReader()->GetEMCALCells(),etaRegion,phiRegion);
870  // Region 0: center of SM ~0.18<|eta|<0.55
871  if ( etaRegion !=0 ) return ;
872 
873  fhEtaPhiFECCorrControl->Fill(fClusterMomentum.Eta(),GetPhi(fClusterMomentum.Phi()));
874 
875  Float_t energy = clus->E();
876  Float_t m02 = clus->GetM02();
877  Float_t m20 = clus->GetM20();
878  Int_t ncells = clus->GetNCells();
879 
880  Int_t absIdCorr[] = {-1,-1,-1,-1};
881  Bool_t correlMaxOK = GetCaloUtils()->GetFECCorrelatedCellAbsId(absIdMax, absIdCorr);
882 
883  //
884  // Correlation to max
885  //
886  Int_t nCorr = 0;
887  Int_t nCellWithWeight = 1;
888  Bool_t near = kFALSE;
889  Bool_t far = kFALSE;
890  // Get second highest energy cell
891  Int_t absId2ndMax = -1;
892  Float_t emax = 0;
893  if(correlMaxOK)
894  {
895  for (Int_t ipos = 0; ipos < ncells; ipos++)
896  {
897  Int_t absId = clus->GetCellsAbsId()[ipos];
898 
899  Float_t eCell = cells->GetCellAmplitude(absId);
900  // consider cells with enough energy weight and not the reference one
901  Float_t weight = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(eCell, energy);
902 
903  if( absId == absIdMax || weight < 0.01 ) continue;
904 
905  nCellWithWeight++;
906 
907  if ( eCell > emax )
908  {
909  emax = eCell;
910  absId2ndMax = absId;
911  }
912 
913  for(Int_t id = 0; id < 4; id++)
914  {
915  if ( absId == absIdCorr[id] )
916  {
917  nCorr++;
918  Int_t diff = TMath::Abs(absId-absIdMax);
919  if ( diff == 1 ) near = kTRUE;
920  else if( diff >= 8 ) far = kTRUE;
921  }
922  }
923  }
924 
925  fhNCellsPerClusterWithWeight->Fill(energy, nCellWithWeight, GetEventWeight());
926  Float_t ratioNcells = nCellWithWeight/(ncells*1.);
927  //printf("E %2.2f, ncells %d, nCellWithWeight %d, ratio %2.2f\n",energy,ncells,nCellWithWeight,ratioNcells);
928  fhNCellsPerClusterRatioWithWeight->Fill(energy, ratioNcells, GetEventWeight());
929 
930  //if ( nCorr || close || far ) printf("nCorr = %d, close %d, far %d; E %2.2f, ncell %d\n",nCorr,close,far,e,nCaloCellsPerCluster);
931 
932  //if ( nCorr > 3) printf("More than 3 correlations: %d\n",nCorr);
933 
934  // Fill histograms assign index
935  Int_t index = -1;
936  if ( nCorr == 0 ) index = 0;
937  else
938  {
939  if ( near && !far ) index = 1;
940  else if (!near && far ) index = 2;
941  else if ( near && far ) index = 3;
942  }
943 
944  if ( index >= 0 )
945  {
946  fhLambda0MaxFECCorrel[index]->Fill(energy, m02, GetEventWeight());
947  fhLambda1MaxFECCorrel[index]->Fill(energy, m20, GetEventWeight());
948 
949  if(nCellWithWeight > 3)
950  {
951  fhLambda0MaxFECCorrelLargeNCells[index]->Fill(energy, m02, GetEventWeight());
952  fhLambda1MaxFECCorrelLargeNCells[index]->Fill(energy, m20, GetEventWeight());
953  }
954  }
955  } // list of correlated towers is found
956  else printf("Max cell correl cells not found\n");
957 
958  //
959  // Any cell correlation, excluding max and their associates
960  //
961  Int_t absIdCorrSecondary [] = {-1,-1,-1,-1};
962  Int_t counter = 0;
963  for (Int_t ipos = 0; ipos < ncells; ipos++)
964  {
965  Int_t absId = clus->GetCellsAbsId()[ipos];
966 
967  // consider cells with enough energy weight and not the reference one or their correlated
968  Float_t weight = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(absId), energy);
969 
970  if( absId == absIdCorr[0] || absId == absIdCorr[1] ||
971  absId == absIdCorr[2] || absId == absIdCorr[3] || weight < 0.01 ) continue;
972 
973  Bool_t correlOK = GetCaloUtils()->GetFECCorrelatedCellAbsId(absId, absIdCorrSecondary);
974 
975  if(!correlOK) continue;
976 
977  for (Int_t ipos2 = ipos+1; ipos2 < ncells; ipos2++)
978  {
979  Int_t absId2 = clus->GetCellsAbsId()[ipos2];
980 
981  // consider cells with enough energy weight and not the reference one or their correlated
982  Float_t weight = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(absId2), energy);
983 
984  if( absId2 == absIdCorr[0] || absId2 == absIdCorr[1] || absId2 == absId ||
985  absId2 == absIdCorr[2] || absId2 == absIdCorr[3] || weight < 0.01 ) continue;
986 
987  for(Int_t id = 0; id < 4; id++)
988  {
989  if ( absId2 == absIdCorrSecondary[id] )
990  counter++;
991  }
992  }
993  }
994 
995  Int_t distCase = 0;
996  if ( counter == 1 ) distCase = 1;
997  else if ( counter >= 2 ) distCase = 2;
998 
999  if ( nCorr > 0 ) distCase+=3;
1000 
1001  fhLambda0FECCorrel[distCase]->Fill(energy, m02, GetEventWeight());
1002  fhLambda1FECCorrel[distCase]->Fill(energy, m20, GetEventWeight());
1003  if(nCellWithWeight > 3)
1004  {
1005  fhLambda0FECCorrelLargeNCells[distCase]->Fill(energy, m02, GetEventWeight());
1006  fhLambda1FECCorrelLargeNCells[distCase]->Fill(energy, m20, GetEventWeight());
1007  }
1008 
1009  //
1010  // Compare pair and odd eta
1011  //
1012  Int_t idMax = -1;
1013  Int_t idMax2 = -1;
1014  Bool_t pairCol = kFALSE;
1015  for(Int_t id = 0; id < 4; id++)
1016  {
1017  if ( absIdMax == absIdCorr[id]) idMax = id;
1018  else if( absId2ndMax == absIdCorr[id]) idMax2 = id;
1019  }
1020 
1021  if(idMax < 0)
1022  {
1023  printf("something wrong, absId max not found in correlated");
1024  return;
1025  }
1026 
1027  if(idMax == 0 || idMax == 3) pairCol = kTRUE; // n and n+8
1028  // else n+1 and n+9
1029 
1030  if(pairCol)
1031  {
1032  fhLambda0MaxFECPair->Fill(energy, m02, GetEventWeight());
1033  fhLambda1MaxFECPair->Fill(energy, m20, GetEventWeight());
1034  if ( idMax2 >=0 )
1035  {
1036  fhLambda0Max2FECPair->Fill(energy, m02, GetEventWeight());
1037  fhLambda1Max2FECPair->Fill(energy, m20, GetEventWeight());
1038  }
1039  }
1040  else
1041  {
1042  fhLambda0MaxFECOdd->Fill(energy, m02, GetEventWeight());
1043  fhLambda1MaxFECOdd->Fill(energy, m20, GetEventWeight());
1044  if ( idMax2 >=0 )
1045  {
1046  fhLambda0Max2FECOdd->Fill(energy, m02, GetEventWeight());
1047  fhLambda1Max2FECOdd->Fill(energy, m20, GetEventWeight());
1048  }
1049  }
1050 
1051  if(nCellWithWeight > 3)
1052  {
1053  if(pairCol)
1054  {
1055  fhLambda0MaxFECPairLargeNCells->Fill(energy, m02, GetEventWeight());
1056  fhLambda1MaxFECPairLargeNCells->Fill(energy, m20, GetEventWeight());
1057  if ( idMax2 >=0 )
1058  {
1059  fhLambda0Max2FECPairLargeNCells->Fill(energy, m02, GetEventWeight());
1060  fhLambda1Max2FECPairLargeNCells->Fill(energy, m20, GetEventWeight());
1061  }
1062  }
1063  else
1064  {
1065  fhLambda0MaxFECOddLargeNCells->Fill(energy, m02, GetEventWeight());
1066  fhLambda1MaxFECOddLargeNCells->Fill(energy, m20, GetEventWeight());
1067  if ( idMax2 >=0 )
1068  {
1069  fhLambda0Max2FECOddLargeNCells->Fill(energy, m02, GetEventWeight());
1070  fhLambda1Max2FECOddLargeNCells->Fill(energy, m20, GetEventWeight());
1071  }
1072  }
1073  }
1074 
1075 }
1076 
1077 //_____________________________________________________________________________________
1078 // Study the shape of the cluster in cell units terms.
1079 //_____________________________________________________________________________________
1080 void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax,
1081  Bool_t goodCluster)
1082 {
1083  // No use to study clusters with less than 4 cells
1084  if( clus->GetNCells() <= 3 ) return;
1085 
1086  Int_t dIeta = 0;
1087  Int_t dIphi = 0;
1088 
1089  Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
1090  Int_t smMax = GetModuleNumberCellIndexes(absIdMax,GetCalorimeter(), ietaMax, iphiMax, rcuMax);
1091 
1092  for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
1093  {
1094  Int_t absId = clus->GetCellsAbsId()[ipos];
1095 
1096  Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
1097  Int_t sm = GetModuleNumberCellIndexes(absId,GetCalorimeter(), ieta, iphi, rcu);
1098 
1099  if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
1100 
1101  if(smMax==sm)
1102  {
1103  if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
1104  }
1105  else
1106  {
1107  Int_t ietaShift = ieta;
1108  Int_t ietaMaxShift = ietaMax;
1109  if (ieta > ietaMax) ietaMaxShift+=48;
1110  else ietaShift +=48;
1111  if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
1112  }
1113  }// Fill cell-cluster histogram loop
1114 
1115  Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
1116 
1117  if(goodCluster)
1118  {
1119  // Was cluster matched?
1120  Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
1121 
1122  if (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta, dIphi, GetEventWeight());
1123  else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta, dIphi, GetEventWeight());
1124  else fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta, dIphi, GetEventWeight());
1125 
1126  fhDeltaIA[matched]->Fill(clus->E(), dIA, GetEventWeight());
1127 
1128  if(clus->E() > 0.5)
1129  {
1130  fhDeltaIAL0 [matched]->Fill(clus->GetM02() , dIA, GetEventWeight());
1131  fhDeltaIAL1 [matched]->Fill(clus->GetM20() , dIA, GetEventWeight());
1132  fhDeltaIANCells[matched]->Fill(clus->GetNCells(), dIA, GetEventWeight());
1133  }
1134 
1135  // Origin of clusters
1136  Int_t nLabel = clus->GetNLabels();
1137  Int_t* labels = clus->GetLabels();
1138 
1139  if(IsDataMC())
1140  {
1141  Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),GetCalorimeter());
1142  if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1143  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) &&
1144  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
1145  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
1146  fhDeltaIAMC[0]->Fill(clus->E(), dIA, GetEventWeight()); // Pure Photon
1147  }
1148  else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1149  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
1150  fhDeltaIAMC[1]->Fill(clus->E(), dIA, GetEventWeight()); // Pure electron
1151  }
1152  else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1153  GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
1154  fhDeltaIAMC[2]->Fill(clus->E(), dIA, GetEventWeight()); // Converted cluster
1155  }
1156  else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
1157  fhDeltaIAMC[3]->Fill(clus->E(), dIA, GetEventWeight()); // Hadrons
1158  }
1159 
1160  } // MC
1161  } // good cluster
1162  else
1163  {
1164  if (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta, dIphi, GetEventWeight());
1165  else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta, dIphi, GetEventWeight());
1166  else fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta, dIphi, GetEventWeight());
1167 
1168  fhBadClusterDeltaIA->Fill(clus->E(), dIA, GetEventWeight());
1169  }
1170 }
1171 
1172 //__________________________________________________________________________________________________________________
1181 //__________________________________________________________________________________________________________________
1182 void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
1183  Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
1184  Double_t tmax)
1185 {
1186  Double_t tof = clus->GetTOF()*1.e9;
1187  if(tof>400) tof-=fConstantTimeShift;
1188 
1189  fhLambda0 ->Fill(clus->E(), clus->GetM02() , GetEventWeight());
1190  fhLambda1 ->Fill(clus->E(), clus->GetM20() , GetEventWeight());
1191  // fhDispersion ->Fill(clus->E(), clus->GetDispersion(), GetEventWeight());
1192 
1193  fhClusterMaxCellDiff ->Fill(clus->E(), maxCellFraction, GetEventWeight());
1194  fhClusterMaxCellECross->Fill(clus->E(), eCrossFrac , GetEventWeight());
1195  fhClusterTimeEnergy ->Fill(clus->E(), tof , GetEventWeight());
1196 
1197  if(fStudyM02Dependence)
1198  {
1199  fhClusterMaxCellDiffM02 ->Fill(clus->E(), maxCellFraction, clus->GetM02(), GetEventWeight());
1200  fhClusterMaxCellECrossM02->Fill(clus->E(), eCrossFrac , clus->GetM02(), GetEventWeight());
1201  fhClusterTimeEnergyM02 ->Fill(clus->E(), tof , clus->GetM02(), GetEventWeight());
1202  }
1203 
1204  if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
1205 
1206  Int_t nModule = GetModuleNumber(clus);
1207  Int_t nCaloCellsPerCluster = clus->GetNCells();
1208 
1209  // Clusters in event time difference
1210  for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
1211  {
1212  AliVCluster* clus2 = (AliVCluster*) caloClusters->At(iclus2);
1213 
1214  if( clus->GetID() == clus2->GetID() ) continue;
1215 
1216  if( clus->GetM02() > 0.01 && clus2->GetM02() > 0.01 )
1217  {
1218  Int_t nModule2 = GetModuleNumber(clus2);
1219 
1220  Double_t tof2 = clus2->GetTOF()*1.e9;
1221  if(tof2>400) tof2-=fConstantTimeShift;
1222 
1223  fhClusterPairDiffTimeE ->Fill(clus->E(), tof-tof2, GetEventWeight());
1224 
1225  if ( nModule2 == nModule )
1226  fhClusterPairDiffTimeESameMod->Fill(clus->E(), tof-tof2, GetEventWeight());
1227  }
1228  }
1229 
1230  if(nCaloCellsPerCluster > 1)
1231  {
1232  // Check time of cells respect to max energy cell
1233 
1234 // if(fFillAllCellTimeHisto)
1235 // {
1236 // // Get some time averages
1237 // Double_t timeAverages[2] = {0.,0.};
1238 // CalculateAverageTime(clus, cells, timeAverages);
1239 //
1240 // fhClusterMaxCellDiffAverageTime ->Fill(clus->E(), tmax-timeAverages[0], GetEventWeight());
1241 // fhClusterMaxCellDiffWeightedTime ->Fill(clus->E(), tmax-timeAverages[1], GetEventWeight());
1242 // }
1243 
1244  for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
1245  {
1246  Int_t absId = clus->GetCellsAbsId()[ipos];
1247  if( absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01 ) continue;
1248 
1249  Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
1250  fhClusterMaxCellCloseCellRatio->Fill(clus->E(), frac, GetEventWeight());
1251  fhClusterMaxCellCloseCellDiff ->Fill(clus->E(), cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId), GetEventWeight());
1252 
1253  if(fStudyM02Dependence)
1254  {
1255  fhClusterMaxCellCloseCellRatioM02->Fill(clus->E(), frac, clus->GetM02(), GetEventWeight());
1256  fhClusterMaxCellCloseCellDiffM02 ->Fill(clus->E(), cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId), clus->GetM02(),GetEventWeight());
1257  }
1258 
1259  if(fFillAllCellTimeHisto)
1260  {
1261  Double_t time = cells->GetCellTime(absId);
1262  GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1263 
1264  Float_t diff = (tmax-(time*1.0e9-fConstantTimeShift));
1265  fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff, GetEventWeight());
1266  if(fStudyM02Dependence)
1267  fhCellTimeSpreadRespectToCellMaxM02->Fill(clus->E(), diff, clus->GetM02(), GetEventWeight());
1268 
1269  if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId, GetEventWeight());
1270  }
1271 
1272  } // Fill cell-cluster histogram loop
1273 
1274  } // Check time and energy of cells respect to max energy cell if cluster of more than 1 cell
1275 
1276  Float_t e = fClusterMomentum.E();
1277  Float_t pt = fClusterMomentum.Pt();
1278  Float_t eta = fClusterMomentum.Eta();
1279  Float_t phi = GetPhi(fClusterMomentum.Phi());
1280 
1281  AliDebug(1,Form("cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f",e,pt,eta,phi*TMath::RadToDeg()));
1282 
1283  fhE ->Fill(e, GetEventWeight());
1284  if(nModule >=0 && nModule < fNModules)
1285  {
1286  fhEMod ->Fill(e, nModule, GetEventWeight());
1287  fhEWeirdMod->Fill(e, nModule, GetEventWeight()); // different binning
1288  }
1289 
1290  fhPt ->Fill(pt , GetEventWeight());
1291  fhPhi ->Fill(phi, GetEventWeight());
1292  fhEta ->Fill(eta, GetEventWeight());
1293 
1294  if(fFillEBinAcceptanceHisto)
1295  {
1296  Int_t icol = -1, irow = -1, iRCU = -1, icolAbs = -1, irowAbs = -1;
1297  GetModuleNumberCellIndexesAbsCaloMap(absIdMax,GetCalorimeter(), icol, irow, iRCU, icolAbs, irowAbs);
1298 
1299  for(Int_t ie = 0; ie < fNEBinCuts; ie++)
1300  {
1301  if( e >= fEBinCuts[ie] && e < fEBinCuts[ie+1] )
1302  {
1303  fhEBinClusterEtaPhi[ie]->Fill(eta,phi,GetEventWeight()) ;
1304 
1305  fhEBinClusterColRow[ie]->Fill(icolAbs,irowAbs,GetEventWeight()) ;
1306  }
1307  }
1308  }
1309  else if ( fFillAllTH3 ) fhEtaPhiE->Fill(eta, phi, e, GetEventWeight());
1310  else if ( e > 0.5 ) fhEtaPhi ->Fill(eta, phi, GetEventWeight());
1311 
1312  // Cells per cluster
1313  fhNCellsPerCluster->Fill(e, nCaloCellsPerCluster, GetEventWeight());
1314 
1315  if(fStudyM02Dependence)
1316  fhNCellsPerClusterM02->Fill(e, nCaloCellsPerCluster, clus->GetM02(),GetEventWeight());
1317 
1318  if(e > 100)
1319  fhNCellsPerClusterWeirdMod->Fill(nCaloCellsPerCluster, nModule, GetEventWeight());
1320 
1321  // Position
1322  if(fFillAllPosHisto2)
1323  {
1324  Float_t pos[3] ;
1325  clus->GetPosition(pos);
1326 
1327  fhXE ->Fill(pos[0], e, GetEventWeight());
1328  fhYE ->Fill(pos[1], e, GetEventWeight());
1329  fhZE ->Fill(pos[2], e, GetEventWeight());
1330 
1331  fhXYZ ->Fill(pos[0], pos[1], pos[2], GetEventWeight());
1332 
1333  fhXNCells->Fill(pos[0], nCaloCellsPerCluster, GetEventWeight());
1334  fhYNCells->Fill(pos[1], nCaloCellsPerCluster, GetEventWeight());
1335  fhZNCells->Fill(pos[2], nCaloCellsPerCluster, GetEventWeight());
1336 
1337  Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
1338 
1339  fhRE ->Fill(rxyz, e , GetEventWeight());
1340  fhRNCells->Fill(rxyz, nCaloCellsPerCluster, GetEventWeight());
1341  }
1342 
1343  if( nModule >= 0 && nModule < fNModules ) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster, GetEventWeight());
1344 }
1345 
1346 //____________________________________________________________________________
1358 //____________________________________________________________________________
1359 void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters,
1360  AliVCaloCells* cells)
1361 {
1362  Int_t nLabel = 0 ;
1363  Int_t *labels = 0x0;
1364  Int_t nCaloClusters = caloClusters->GetEntriesFast() ;
1365  Int_t nCaloClustersAccepted = 0 ;
1366  Int_t nCaloCellsPerCluster = 0 ;
1367  Bool_t matched = kFALSE;
1368  Int_t nModule =-1 ;
1369 
1370  // Get vertex for photon momentum calculation and event selection
1371  Double_t v[3] = {0,0,0}; //vertex ;
1372 //GetReader()->GetVertex(v);
1373 
1374  Int_t *nClustersInModule = new Int_t [fNModules];
1375  Float_t *energyInModule = new Float_t[fNModules];
1376  for(Int_t imod = 0; imod < fNModules; imod++ )
1377  {
1378  nClustersInModule[imod] = 0;
1379  energyInModule [imod] = 0;
1380  }
1381 
1382  AliDebug(1,Form("In %s there are %d clusters", GetCalorimeterString().Data(), nCaloClusters));
1383 
1384  // Loop over CaloClusters
1385  for(Int_t iclus = 0; iclus < nCaloClusters; iclus++)
1386  {
1387  AliDebug(1,Form("Cluster: %d/%d, data %d",iclus+1,nCaloClusters,GetReader()->GetDataType()));
1388 
1389  AliVCluster* clus = (AliVCluster*) caloClusters->At(iclus);
1390 
1391  // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
1392  Float_t maxCellFraction = 0.;
1393  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
1394 
1395  // Cut on time of clusters
1396  Double_t tof = clus->GetTOF()*1.e9;
1397  if(tof>400) tof-=fConstantTimeShift;
1398 
1399  if( tof < fTimeCutMin || tof > fTimeCutMax )
1400  {
1401  AliDebug(1,Form("Remove cluster with TOF %f",tof));
1402  continue;
1403  }
1404 
1405  // Get cluster kinematics
1406  clus->GetMomentum(fClusterMomentum,v);
1407 
1408  // Check only certain regions
1409  Bool_t in = kTRUE;
1410  if(IsFiducialCutOn()) in = GetFiducialCut()->IsInFiducialCut(fClusterMomentum.Eta(),fClusterMomentum.Phi(),GetCalorimeter()) ;
1411  if(!in) continue;
1412 
1413  // MC labels
1414  nLabel = clus->GetNLabels();
1415  labels = clus->GetLabels();
1416 
1417  // SuperModule number of cluster
1418  nModule = GetModuleNumber(clus);
1419 
1420  // Cells per cluster
1421  nCaloCellsPerCluster = clus->GetNCells();
1422 
1423  // Cluster mathed with track?
1424  matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
1425 
1426  // Get time of max cell
1427  Double_t tmax = cells->GetCellTime(absIdMax);
1428  GetCaloUtils()->RecalibrateCellTime(tmax, GetCalorimeter(), absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1429  tmax*=1.e9;
1430  tmax-=fConstantTimeShift;
1431  //
1432  // Fill histograms related to single cluster
1433  //
1434 
1435  // Fill some histograms before applying the exotic cell / bad map cut
1436  if(fStudyBadClusters)
1437  {
1438  fhNCellsPerClusterNoCut->Fill(clus->E(), nCaloCellsPerCluster, GetEventWeight());
1439 
1440  if(nModule >=0 && nModule < fNModules)
1441  {
1442  fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster, GetEventWeight());
1443  if(clus->E() > 100) fhNCellsPerClusterWeirdModNoCut->Fill(nCaloCellsPerCluster, nModule, GetEventWeight());
1444  }
1445 
1446  fhClusterMaxCellDiffNoCut->Fill(clus->E(), maxCellFraction, GetEventWeight());
1447  }
1448 
1449  Float_t ampMax = cells->GetCellAmplitude(absIdMax);
1450  GetCaloUtils()->RecalibrateCellAmplitude(ampMax,GetCalorimeter(), absIdMax);
1451 
1452  if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
1453 
1454  // Check bad clusters if requested and rejection was not on
1455  Bool_t goodCluster = IsGoodCluster(absIdMax, clus->GetM02(), nCaloCellsPerCluster, cells);
1456 
1457  Float_t eCrossFrac = 0;
1458  if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
1459 
1460  AliDebug(1,Form("Accept cluster? %d",goodCluster));
1461 
1462  if(!goodCluster)
1463  {
1464  if ( fStudyBadClusters ) BadClusterHistograms(clus, caloClusters,
1465  cells, absIdMax,
1466  maxCellFraction,
1467  eCrossFrac, tmax);
1468  continue;
1469  }
1470 
1471  //
1472  ClusterHistograms(clus, caloClusters, cells, absIdMax,
1473  maxCellFraction, eCrossFrac, tmax);
1474 
1475  nCaloClustersAccepted++;
1476 
1477  //
1478  if(fStudyFECCorrelation) ChannelCorrelationInFEC(clus, cells, matched, absIdMax);
1479 
1480  //
1481  nModule = GetModuleNumber(clus);
1482  if(nModule >=0 && nModule < fNModules && fClusterMomentum.E() > 2*fCellAmpMin)
1483  {
1484  nClustersInModule[nModule]++;
1485  if(clus->E() > 0.5)
1486  energyInModule [nModule] += clus->E();
1487  }
1488 
1489  // Cluster weights
1490  if(fStudyWeight) WeightHistograms(clus, cells);
1491 
1492  // Cells in cluster position
1493  if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
1494 
1495  // Fill histograms related to single cluster, mc vs data
1496  Int_t mcOK = kFALSE;
1497  Int_t pdg = -1;
1498  if(IsDataMC() && nLabel > 0 && labels)
1499  mcOK = ClusterMCHistograms(matched, labels, nLabel, pdg);
1500 
1501  // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
1502  if( matched && fFillAllTMHisto)
1503  ClusterMatchedWithTrackHistograms(clus,mcOK,pdg);
1504 
1505  // Invariant mass
1506  // Try to reduce background with a mild shower shape cut and no more than 1 maxima
1507  // in cluster and remove low energy clusters
1508 
1509  if ( fFillAllPi0Histo
1510  && nCaloClusters > 1
1511  && nCaloCellsPerCluster > 1
1512  && GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1
1513  && clus->GetM02() < fInvMassMaxM02Cut
1514  && clus->GetM02() > fInvMassMinM02Cut
1515  && clus->E() > fInvMassMinECut
1516  && clus->E() < fInvMassMaxECut
1517  )
1518  InvariantMassHistograms(iclus, nModule, caloClusters,cells);
1519 
1520  } // Cluster loop
1521 
1522  // Number of clusters histograms
1523  if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted, GetEventWeight());
1524 
1525  // Number of clusters per module
1526  for(Int_t imod = 0; imod < fNModules; imod++ )
1527  {
1528  AliDebug(1,Form("Module %d calo %s clusters %d, sum E %f", imod, GetCalorimeterString().Data(), nClustersInModule[imod], energyInModule[imod]));
1529 
1530  fhNClustersMod ->Fill(nClustersInModule[imod], imod, GetEventWeight());
1531  fhSumClustersEnergyMod->Fill(energyInModule [imod], imod, GetEventWeight());
1532 
1533  fhNClustersSumEnergyPerMod[imod]->Fill(energyInModule[imod], nClustersInModule[imod], GetEventWeight());
1534  }
1535 
1536  delete [] nClustersInModule;
1537  delete [] energyInModule;
1538 }
1539 
1540 //__________________________________________________________________________________
1543 //__________________________________________________________________________________
1544 Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(Bool_t matched,const Int_t * labels,
1545  Int_t nLabels, Int_t & pdg )
1546 {
1547  if(!labels || nLabels<=0)
1548  {
1549  AliWarning(Form("Strange, labels array %p, n labels %d", labels,nLabels));
1550  return kFALSE;
1551  }
1552 
1553  AliDebug(1,Form("Primaries: nlabels %d",nLabels));
1554 
1555  Float_t e = fClusterMomentum.E();
1556  Float_t eta = fClusterMomentum.Eta();
1557  Float_t phi = fClusterMomentum.Phi();
1558  if(phi < 0) phi +=TMath::TwoPi();
1559 
1560  AliAODMCParticle * aodprimary = 0x0;
1561  TParticle * primary = 0x0;
1562 
1563  // Play with the MC stack if available
1564  Int_t label = labels[0];
1565 
1566  if(label < 0)
1567  {
1568  AliDebug(1,Form(" *** bad label ***: label %d", label));
1569  return kFALSE;
1570  }
1571 
1572  Int_t pdg0 =-1; Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1573  Float_t vxMC = 0; Float_t vyMC = 0;
1574  Float_t eMC = 0; //Float_t ptMC= 0;
1575  Float_t phiMC = 0; Float_t etaMC = 0;
1576  Int_t charge = 0;
1577 
1578  // Check the origin.
1579  Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),GetCalorimeter());
1580 
1581  if ( GetReader()->ReadStack() &&
1582  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1583  { // If MC stack and known tag
1584 
1585  if( label >= GetMCStack()->GetNtrack())
1586  {
1587  AliDebug(1,Form("*** large label ***: label %d, n tracks %d", label, GetMCStack()->GetNtrack()));
1588  return kFALSE;
1589  }
1590 
1591  primary = GetMCStack()->Particle(label);
1592  iMother = label;
1593  pdg0 = TMath::Abs(primary->GetPdgCode());
1594  pdg = pdg0;
1595  status = primary->GetStatusCode();
1596  vxMC = primary->Vx();
1597  vyMC = primary->Vy();
1598  iParent = primary->GetFirstMother();
1599 
1600  AliDebug(1,"Cluster most contributing mother:");
1601  AliDebug(1,Form("\t Mother label %d, pdg %d, %s, status %d, parent %d",iMother, pdg0, primary->GetName(),status, iParent));
1602 
1603 
1604  // Get final particle, no conversion products
1605  if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1606  {
1607  // Get the parent
1608  primary = GetMCStack()->Particle(iParent);
1609  pdg = TMath::Abs(primary->GetPdgCode());
1610 
1611  AliDebug(2,"Converted cluster!. Find before conversion:");
1612 
1613  while((pdg == 22 || pdg == 11) && status != 1)
1614  {
1615  Int_t iMotherOrg = iMother;
1616  iMother = iParent;
1617  primary = GetMCStack()->Particle(iMother);
1618  status = primary->GetStatusCode();
1619  pdg = TMath::Abs(primary->GetPdgCode());
1620  iParent = primary->GetFirstMother();
1621 
1622  // If gone too back and non stable, assign the decay photon/electron
1623  // there are other possible decays, ignore them for the moment
1624  if(pdg==111 || pdg==221)
1625  {
1626  primary = GetMCStack()->Particle(iMotherOrg);
1627  break;
1628  }
1629 
1630  if( iParent < 0 )
1631  {
1632  iParent = iMother;
1633  break;
1634  }
1635 
1636  AliDebug(2,Form("\t pdg %d, index %d, %s, status %d",pdg, iMother, primary->GetName(),status));
1637  }
1638 
1639  AliDebug(1,"Converted Cluster mother before conversion:");
1640  AliDebug(1,Form("\t Mother label %d, pdg %d, %s, status %d, parent %d",iMother, pdg, primary->GetName(), status, iParent));
1641 
1642  }
1643 
1644  // Overlapped pi0 (or eta, there will be very few), get the meson
1645  if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1646  GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1647  {
1648  AliDebug(2,"Overlapped Meson decay!, Find it:");
1649 
1650  while(pdg != 111 && pdg != 221)
1651  {
1652  //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1653  iMother = iParent;
1654  primary = GetMCStack()->Particle(iMother);
1655  status = primary->GetStatusCode();
1656  pdg = TMath::Abs(primary->GetPdgCode());
1657  iParent = primary->GetFirstMother();
1658 
1659  if( iParent < 0 ) break;
1660 
1661  AliDebug(2,Form("\t pdg %d, %s, index %d",pdg, primary->GetName(),iMother));
1662 
1663  if(iMother==-1)
1664  {
1665  AliWarning("Tagged as Overlapped photon but meson not found, why?");
1666  //break;
1667  }
1668  }
1669 
1670  AliDebug(2,Form("Overlapped %s decay, label %d",primary->GetName(),iMother));
1671  }
1672 
1673  eMC = primary->Energy();
1674  //ptMC = primary->Pt();
1675  phiMC = primary->Phi();
1676  etaMC = primary->Eta();
1677  pdg = TMath::Abs(primary->GetPdgCode());
1678  charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1679 
1680  }
1681  else if( GetReader()->ReadAODMCParticles() &&
1682  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1683  {
1684  // If MC AOD and known tag
1685  // Get the list of MC particles
1686  if(!GetReader()->GetAODMCParticles())
1687  AliFatal("MCParticles not available!");
1688 
1689  aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(label);
1690  iMother = label;
1691  pdg0 = TMath::Abs(aodprimary->GetPdgCode());
1692  pdg = pdg0;
1693  status = aodprimary->IsPrimary();
1694  vxMC = aodprimary->Xv();
1695  vyMC = aodprimary->Yv();
1696  iParent = aodprimary->GetMother();
1697 
1698  AliDebug(1,"Cluster most contributing mother:");
1699  AliDebug(1,Form("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d",
1700  iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent));
1701 
1702  //Get final particle, no conversion products
1703  if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) )
1704  {
1705  AliDebug(2,"Converted cluster!. Find before conversion:");
1706 
1707  // Get the parent
1708  aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iParent);
1709  pdg = TMath::Abs(aodprimary->GetPdgCode());
1710 
1711  while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1712  {
1713  Int_t iMotherOrg = iMother;
1714  iMother = iParent;
1715  aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1716  status = aodprimary->IsPrimary();
1717  iParent = aodprimary->GetMother();
1718  pdg = TMath::Abs(aodprimary->GetPdgCode());
1719 
1720  // If gone too back and non stable, assign the decay photon/electron
1721  // there are other possible decays, ignore them for the moment
1722  if( pdg == 111 || pdg == 221 )
1723  {
1724  aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMotherOrg);
1725  break;
1726  }
1727 
1728  if( iParent < 0 )
1729  {
1730  iParent = iMother;
1731  break;
1732  }
1733 
1734  AliDebug(2,Form("\t pdg %d, index %d, Primary? %d, Physical Primary? %d",
1735  pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()));
1736  }
1737 
1738  AliDebug(1,"Converted Cluster mother before conversion:");
1739  AliDebug(1,Form("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d",
1740  iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()));
1741 
1742  }
1743 
1744  //Overlapped pi0 (or eta, there will be very few), get the meson
1745  if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1746  GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1747  {
1748  AliDebug(2,Form("Overlapped Meson decay!, Find it: PDG %d, mom %d",pdg, iMother));
1749 
1750  while(pdg != 111 && pdg != 221)
1751  {
1752  iMother = iParent;
1753  aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1754  status = aodprimary->IsPrimary();
1755  iParent = aodprimary->GetMother();
1756  pdg = TMath::Abs(aodprimary->GetPdgCode());
1757 
1758  if( iParent < 0 ) break;
1759 
1760  AliDebug(2,Form("\t pdg %d, index %d",pdg, iMother));
1761 
1762  if(iMother==-1)
1763  {
1764  AliWarning("Tagged as Overlapped photon but meson not found, why?");
1765  //break;
1766  }
1767  }
1768 
1769  AliDebug(2,Form("Overlapped %s decay, label %d",aodprimary->GetName(),iMother));
1770  }
1771 
1772  status = aodprimary->IsPrimary();
1773  eMC = aodprimary->E();
1774  //ptMC = aodprimary->Pt();
1775  phiMC = aodprimary->Phi();
1776  etaMC = aodprimary->Eta();
1777  pdg = TMath::Abs(aodprimary->GetPdgCode());
1778  charge = aodprimary->Charge();
1779  }
1780 
1781  //Float_t vz = primary->Vz();
1782  Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1783  if( ( pdg == 22 || TMath::Abs(pdg) == 11 ) && status != 1 )
1784  {
1785  fhEMVxyz ->Fill(vxMC, vyMC, GetEventWeight());//,vz);
1786  fhEMR ->Fill(e , rVMC, GetEventWeight());
1787  }
1788 
1789  //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1790  //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1791  //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1792 
1793  // Overlapped pi0 (or eta, there will be very few)
1794  Int_t mcIndex = -1;
1795  if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0 ) )
1796  {
1797  mcIndex = kmcPi0;
1798  } // Overlapped pizero decay
1799  else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta ) )
1800  {
1801  mcIndex = kmcEta;
1802  } // Overlapped eta decay
1803  else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton ) )
1804  {
1805  if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1806  mcIndex = kmcPhotonConv ;
1807  else
1808  mcIndex = kmcPhoton ;
1809  } // photon
1810  else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) )
1811  {
1812  mcIndex = kmcElectron;
1813  fhEMVxyz->Fill(vxMC, vyMC, GetEventWeight());//,vz);
1814  fhEMR ->Fill(e , rVMC, GetEventWeight());
1815  }
1816  else if(charge == 0)
1817  {
1818  mcIndex = kmcNeHadron;
1819  fhHaVxyz->Fill(vxMC, vyMC, GetEventWeight());//,vz);
1820  fhHaR ->Fill(e , rVMC, GetEventWeight());
1821  }
1822  else if(charge!=0)
1823  {
1824  mcIndex = kmcChHadron;
1825  fhHaVxyz->Fill(vxMC, vyMC, GetEventWeight());//,vz);
1826  fhHaR ->Fill(e , rVMC, GetEventWeight());
1827  }
1828 
1829  //printf("mc index %d\n",mcIndex);
1830 
1831  if( mcIndex >= 0 && mcIndex < 7 && e > 0.5 && eMC > 0.5)
1832  {
1833  fhRecoMCE [mcIndex][(matched)]->Fill(e , eMC , GetEventWeight());
1834  fhRecoMCEta [mcIndex][(matched)]->Fill(eta, etaMC , GetEventWeight());
1835  fhRecoMCPhi [mcIndex][(matched)]->Fill(phi, phiMC , GetEventWeight());
1836  fhRecoMCRatioE [mcIndex][(matched)]->Fill(e , e/eMC , GetEventWeight());
1837  fhRecoMCDeltaE [mcIndex][(matched)]->Fill(e , eMC-e , GetEventWeight());
1838  fhRecoMCDeltaPhi[mcIndex][(matched)]->Fill(e , phiMC-phi, GetEventWeight());
1839  fhRecoMCDeltaEta[mcIndex][(matched)]->Fill(e , etaMC-eta, GetEventWeight());
1840  }
1841 
1842  if( primary || aodprimary ) return kTRUE ;
1843  else return kFALSE;
1844 }
1845 
1846 //_________________________________________________________________________________________________________
1851 //_________________________________________________________________________________________________________
1852 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, Bool_t okPrimary, Int_t pdg)
1853 {
1854  Float_t e = fClusterMomentum.E();
1855  Float_t pt = fClusterMomentum.Pt();
1856  Float_t eta = fClusterMomentum.Eta();
1857  Float_t phi = fClusterMomentum.Phi();
1858  if(phi < 0) phi +=TMath::TwoPi();
1859 
1860  fhECharged ->Fill(e , GetEventWeight());
1861  fhPtCharged ->Fill(pt , GetEventWeight());
1862  fhPhiCharged ->Fill(phi, GetEventWeight());
1863  fhEtaCharged ->Fill(eta, GetEventWeight());
1864 
1865  if ( fFillAllTH3 ) fhEtaPhiECharged->Fill(eta, phi, e, GetEventWeight());
1866  else if ( e > 0.5 ) fhEtaPhiCharged ->Fill(eta, phi, GetEventWeight());
1867 
1868  // Study the track and matched cluster if track exists.
1869 
1870  AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1871 
1872  if(!track) return ;
1873 
1874  Double_t tpt = track->Pt();
1875  Double_t tmom = track->P();
1876  Double_t dedx = track->GetTPCsignal();
1877  Int_t nITS = track->GetNcls(0);
1878  Int_t nTPC = track->GetNcls(1);
1879  Bool_t positive = kFALSE;
1880  if(track) positive = (track->Charge()>0);
1881 
1882  // Residuals
1883  Float_t deta = clus->GetTrackDz();
1884  Float_t dphi = clus->GetTrackDx();
1885  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1886  Int_t nModule = GetModuleNumber(clus);
1887 
1888  if( TMath::Abs(dphi) < 999 )
1889  {
1890  if(positive)
1891  {
1892  fhTrackMatchedDEtaPos->Fill(e, deta, GetEventWeight());
1893  fhTrackMatchedDPhiPos->Fill(e, dphi, GetEventWeight());
1894 
1895  if(e > 0.5)
1896  {
1897  fhTrackMatchedDEtaPosMod ->Fill(deta, nModule, GetEventWeight());
1898  fhTrackMatchedDPhiPosMod ->Fill(dphi, nModule, GetEventWeight());
1899  fhTrackMatchedDEtaDPhiPos->Fill(deta, dphi , GetEventWeight());
1900  }
1901  }
1902  else
1903  {
1904  fhTrackMatchedDEtaNeg->Fill(e, deta, GetEventWeight());
1905  fhTrackMatchedDPhiNeg->Fill(e, dphi, GetEventWeight());
1906 
1907  if(e > 0.5)
1908  {
1909  fhTrackMatchedDEtaNegMod ->Fill(deta, nModule, GetEventWeight());
1910  fhTrackMatchedDPhiNegMod ->Fill(dphi, nModule, GetEventWeight());
1911  fhTrackMatchedDEtaDPhiNeg->Fill(deta, dphi , GetEventWeight());
1912  }
1913  }
1914  }
1915 
1916  Double_t eOverP = e/tmom;
1917  fh1EOverP->Fill(tpt, eOverP, GetEventWeight());
1918  if(e > 0.5 && tpt > 0.5) fh1EOverPMod->Fill(eOverP, nModule, GetEventWeight());
1919 
1920  if(dR < 0.02)
1921  {
1922  fh1EOverPR02->Fill(tpt, eOverP, GetEventWeight());
1923  if(e > 0.5 && tpt > 0.5) fh1EOverPR02Mod->Fill(eOverP, nModule, GetEventWeight());
1924 
1925  if(dedx > 60 && dedx < 100)
1926  {
1927  fh1EleEOverP->Fill(tpt, eOverP, GetEventWeight());
1928  if(e > 0.5 && tpt > 0.5) fh1EleEOverPMod->Fill(eOverP, nModule, GetEventWeight());
1929  }
1930  }
1931 
1932  fh2dR->Fill(e, dR, GetEventWeight());
1933  fh2MatchdEdx->Fill(tmom, dedx, GetEventWeight());
1934 
1935  if(e > 0.5 && tmom > 0.5)
1936  {
1937  fh2dRMod->Fill(dR, nModule, GetEventWeight());
1938  fh2MatchdEdxMod->Fill(dedx, nModule, GetEventWeight());
1939  }
1940 
1941  if(dR < 0.02 && eOverP > 0.6 && eOverP < 1.2
1942  && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1943  {
1944  fh2EledEdx->Fill(tmom, dedx, GetEventWeight());
1945  if(e > 0.5 && tmom > 0.5) fh2EledEdxMod->Fill(dedx, nModule, GetEventWeight());
1946  }
1947 
1948  if(IsDataMC() && okPrimary)
1949  {
1950  Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1951 
1952  if(TMath::Abs(pdg) == 11)
1953  {
1954  fhMCEle1EOverP ->Fill(tpt , eOverP, GetEventWeight());
1955  fhMCEle1dR ->Fill(dR , GetEventWeight());
1956  fhMCEle2MatchdEdx->Fill(tmom, dedx , GetEventWeight());
1957 
1958  if(dR < 0.02)
1959  {
1960  fhMCEle1EOverPR02->Fill(tpt, eOverP, GetEventWeight());
1961  if(dedx > 60 && dedx < 100) fhMCEle1EleEOverP->Fill(tpt, eOverP, GetEventWeight());
1962  }
1963  }
1964  else if(charge!=0)
1965  {
1966  fhMCChHad1EOverP ->Fill(tpt , eOverP, GetEventWeight());
1967  fhMCChHad1dR ->Fill(dR , GetEventWeight());
1968  fhMCChHad2MatchdEdx->Fill(tmom, dedx , GetEventWeight());
1969 
1970  if(dR < 0.02)
1971  {
1972  fhMCChHad1EOverPR02->Fill(tpt, eOverP, GetEventWeight());
1973  if(dedx > 60 && dedx < 100) fhMCChHad1EleEOverP->Fill(tpt, eOverP, GetEventWeight());
1974  }
1975  }
1976  else if(charge == 0)
1977  {
1978  fhMCNeutral1EOverP ->Fill(tpt , eOverP, GetEventWeight());
1979  fhMCNeutral1dR ->Fill(dR , GetEventWeight());
1980  fhMCNeutral2MatchdEdx->Fill(tmom, dedx , GetEventWeight());
1981 
1982  if(dR < 0.02)
1983  {
1984  fhMCNeutral1EOverPR02->Fill(tpt, eOverP, GetEventWeight());
1985  if(dedx > 60 && dedx < 100) fhMCNeutral1EleEOverP->Fill(tpt, eOverP, GetEventWeight());
1986  }
1987  }
1988  } // DataMC
1989 }
1990 
1991 //___________________________________
1994 //___________________________________
1995 void AliAnaCalorimeterQA::Correlate()
1996 {
1997  // Clusters arrays
1998  TObjArray * caloClustersEMCAL = GetEMCALClusters();
1999  TObjArray * caloClustersPHOS = GetPHOSClusters();
2000 
2001  if(!caloClustersEMCAL || !caloClustersPHOS)
2002  {
2003  AliDebug(1,Form("PHOS (%p) or EMCAL (%p) clusters array not available, do not correlate",caloClustersPHOS,caloClustersEMCAL));
2004  return ;
2005  }
2006 
2007  // Cells arrays
2008  AliVCaloCells * cellsEMCAL = GetEMCALCells();
2009  AliVCaloCells * cellsPHOS = GetPHOSCells();
2010 
2011  if(!cellsEMCAL || !cellsPHOS)
2012  {
2013  AliDebug(1,Form("PHOS (%p) or EMCAL (%p) cells array ot available, do not correlate",cellsPHOS,cellsEMCAL));
2014  return ;
2015  }
2016 
2017  // Clusters parameters
2018  Int_t nclEMCAL = 0;
2019  Int_t nclDCAL = 0;
2020  Int_t nclPHOS = 0;
2021 
2022  Float_t sumClusterEnergyEMCAL = 0;
2023  Float_t sumClusterEnergyDCAL = 0;
2024  Float_t sumClusterEnergyPHOS = 0;
2025 
2026  Int_t iclus = 0;
2027  Float_t energy = 0;
2028  AliVCluster* cluster = 0;
2029  for(iclus = 0 ; iclus < caloClustersEMCAL->GetEntriesFast() ; iclus++)
2030  {
2031  cluster = (AliVCluster*)caloClustersEMCAL->At(iclus);
2032  Float_t energy = cluster->E();
2033 
2034  if( energy < 0.5 ) continue;
2035 
2036  if(cluster->GetCellsAbsId()[0] < 12288)
2037  {
2038  nclEMCAL++;
2039  sumClusterEnergyEMCAL += energy;
2040  }
2041  else
2042  {
2043  nclDCAL++;
2044  sumClusterEnergyDCAL += energy;
2045  }
2046  }
2047 
2048  for(iclus = 0 ; iclus < caloClustersPHOS ->GetEntriesFast(); iclus++)
2049  {
2050  cluster = (AliVCluster*) caloClustersPHOS->At(iclus);
2051 
2052  energy = cluster->E();
2053 
2054  if( energy < 0.5 ) continue;
2055 
2056  nclPHOS++;
2057  sumClusterEnergyPHOS += energy;
2058  }
2059 
2060  // Cells parameters
2061  Int_t ncellsEMCAL = 0 ;
2062  Int_t ncellsDCAL = 0 ;
2063  Int_t ncellsPHOS = 0;
2064 
2065  Float_t sumCellEnergyEMCAL = 0;
2066  Float_t sumCellEnergyDCAL = 0;
2067  Float_t sumCellEnergyPHOS = 0;
2068  Int_t icell = 0;
2069  for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells() ; icell++)
2070  {
2071  Float_t amp = cellsEMCAL->GetAmplitude(icell);
2072  Int_t cellId = cellsEMCAL->GetCellNumber(icell);
2073 
2074  if (amp < fEMCALCellAmpMin) continue;
2075 
2076  if( cellId < 12288)
2077  {
2078  ncellsEMCAL++;
2079  sumCellEnergyEMCAL += amp;
2080  }
2081  else
2082  {
2083  ncellsDCAL++;
2084  sumCellEnergyDCAL += amp;
2085  }
2086  }
2087 
2088  for(icell = 0 ; icell < cellsPHOS->GetNumberOfCells(); icell++)
2089  {
2090  Float_t amp = cellsPHOS->GetAmplitude(icell);
2091  Int_t cellId = cellsPHOS->GetCellNumber(icell);
2092 
2093  if ( cellId < 0 ) continue ; // CPV
2094 
2095  if ( amp < fPHOSCellAmpMin ) continue;
2096 
2097  ncellsPHOS++;
2098  sumCellEnergyPHOS += amp;
2099  }
2100 
2101  // Fill Histograms
2102 
2103  fhEMCALPHOSCorrNClusters->Fill(nclEMCAL , nclPHOS , GetEventWeight());
2104  fhEMCALPHOSCorrEClusters->Fill(sumClusterEnergyEMCAL, sumClusterEnergyPHOS, GetEventWeight());
2105  fhEMCALPHOSCorrNCells ->Fill(ncellsEMCAL , ncellsPHOS , GetEventWeight());
2106  fhEMCALPHOSCorrECells ->Fill(sumCellEnergyEMCAL , sumCellEnergyPHOS , GetEventWeight());
2107 
2108  fhEMCALDCALCorrNClusters->Fill(nclEMCAL , nclDCAL , GetEventWeight());
2109  fhEMCALDCALCorrEClusters->Fill(sumClusterEnergyEMCAL, sumClusterEnergyDCAL, GetEventWeight());
2110  fhEMCALDCALCorrNCells ->Fill(ncellsEMCAL , ncellsDCAL , GetEventWeight());
2111  fhEMCALDCALCorrECells ->Fill(sumCellEnergyEMCAL , sumCellEnergyDCAL , GetEventWeight());
2112 
2113  fhDCALPHOSCorrNClusters ->Fill(nclDCAL , nclPHOS , GetEventWeight());
2114  fhDCALPHOSCorrEClusters ->Fill(sumClusterEnergyDCAL , sumClusterEnergyPHOS, GetEventWeight());
2115  fhDCALPHOSCorrNCells ->Fill(ncellsDCAL , ncellsPHOS , GetEventWeight());
2116  fhDCALPHOSCorrECells ->Fill(sumCellEnergyDCAL , sumCellEnergyPHOS , GetEventWeight());
2117 
2118  Int_t v0S = GetV0Signal(0) + GetV0Signal(1);
2119  Int_t v0M = GetV0Multiplicity(0) + GetV0Multiplicity(1);
2120  Int_t trM = GetTrackMultiplicity();
2121  Float_t cen = GetEventCentrality();
2122  Float_t ep = GetEventPlaneAngle();
2123 
2124  Int_t ncl = nclPHOS;
2125  Float_t sumClusterEnergy = sumClusterEnergyPHOS;
2126  Int_t ncells = ncellsPHOS;
2127  Float_t sumCellEnergy = sumCellEnergyPHOS;
2128 
2129  if ( GetCalorimeter() == kEMCAL )
2130  {
2131  ncl = nclEMCAL + nclDCAL;
2132  sumClusterEnergy = sumClusterEnergyEMCAL + sumClusterEnergyDCAL;
2133  ncells = ncellsEMCAL + ncellsDCAL;
2134  sumCellEnergy = sumCellEnergyEMCAL + sumCellEnergyDCAL;
2135  }
2136 
2137  fhCaloV0MCorrNClusters ->Fill(v0M, ncl , GetEventWeight());
2138  fhCaloV0MCorrEClusters ->Fill(v0M, sumClusterEnergy, GetEventWeight());
2139  fhCaloV0MCorrNCells ->Fill(v0M, ncells , GetEventWeight());
2140  fhCaloV0MCorrECells ->Fill(v0M, sumCellEnergy , GetEventWeight());
2141 
2142  fhCaloV0SCorrNClusters ->Fill(v0S, ncl , GetEventWeight());
2143  fhCaloV0SCorrEClusters ->Fill(v0S, sumClusterEnergy, GetEventWeight());
2144  fhCaloV0SCorrNCells ->Fill(v0S, ncells , GetEventWeight());
2145  fhCaloV0SCorrECells ->Fill(v0S, sumCellEnergy , GetEventWeight());
2146 
2147  fhCaloTrackMCorrNClusters->Fill(trM, ncl , GetEventWeight());
2148  fhCaloTrackMCorrEClusters->Fill(trM, sumClusterEnergy, GetEventWeight());
2149  fhCaloTrackMCorrNCells ->Fill(trM, ncells , GetEventWeight());
2150  fhCaloTrackMCorrECells ->Fill(trM, sumCellEnergy , GetEventWeight());
2151 
2152  fhCaloCenNClusters ->Fill(cen, ncl , GetEventWeight());
2153  fhCaloCenEClusters ->Fill(cen, sumClusterEnergy, GetEventWeight());
2154  fhCaloCenNCells ->Fill(cen, ncells , GetEventWeight());
2155  fhCaloCenECells ->Fill(cen, sumCellEnergy , GetEventWeight());
2156 
2157  fhCaloEvPNClusters ->Fill(ep , ncl , GetEventWeight());
2158  fhCaloEvPEClusters ->Fill(ep , sumClusterEnergy, GetEventWeight());
2159  fhCaloEvPNCells ->Fill(ep , ncells , GetEventWeight());
2160  fhCaloEvPECells ->Fill(ep , sumCellEnergy , GetEventWeight());
2161 
2162  AliDebug(1,"Correlate():");
2163  AliDebug(1,Form("\t EMCAL: N cells %d, N clusters %d, summed E cells %f, summed E clusters %f",
2164  ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL));
2165  AliDebug(1,Form("\t DCAL : N cells %d, N clusters %d, summed E cells %f, summed E clusters %f",
2166  ncellsDCAL,nclDCAL, sumCellEnergyDCAL,sumClusterEnergyDCAL));
2167  AliDebug(1,Form("\t PHOS : N cells %d, N clusters %d, summed E cells %f, summed E clusters %f",
2168  ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS));
2169  AliDebug(1,Form("\t V0 : Signal %d, Multiplicity %d, Track Multiplicity %d", v0S,v0M,trM));
2170  AliDebug(1,Form("\t centrality : %f, Event plane angle %f", cen,ep));
2171 }
2172 
2173 //_________________________________________________
2175 //_________________________________________________
2176 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
2177 {
2178  TString parList ; //this will be list of parameters used for this analysis.
2179  const Int_t buffersize = 255;
2180  char onePar[buffersize] ;
2181 
2182  snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---:") ;
2183  parList+=onePar ;
2184  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
2185  parList+=onePar ;
2186  snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns;",fTimeCutMin, fTimeCutMax) ;
2187  parList+=onePar ;
2188  snprintf(onePar,buffersize,"Cell Amplitude: PHOS > %2.2f GeV, EMCAL > %2.2f GeV;",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
2189  parList+=onePar ;
2190  snprintf(onePar,buffersize,"N cells per cluster: PHOS >= %d, EMCAL >= %d;",fPHOSClusterNCellMin, fEMCALClusterNCellMin) ;
2191  parList+=onePar ;
2192  snprintf(onePar,buffersize,"Inv. Mass %2.1f < E_cl < %2.1f GeV;",fInvMassMinECut, fInvMassMaxECut) ;
2193  parList+=onePar ;
2194  snprintf(onePar,buffersize,"Inv. Mass %2.1f < M02 < %2.1f GeV;",fInvMassMinM02Cut, fInvMassMaxM02Cut) ;
2195  parList+=onePar ;
2196  snprintf(onePar,buffersize,"Cluster pair opening angle < %2.1f rad;",fInvMassMaxOpenAngle) ;
2197  parList+=onePar ;
2198  snprintf(onePar,buffersize,"Cluster pair time difference < %2.1f rad;",fInvMassMaxTimeDifference) ;
2199  parList+=onePar ;
2200 
2201  //Get parameters set in base class.
2202  //parList += GetBaseParametersList() ;
2203 
2204  //Get parameters set in FiducialCut class (not available yet)
2205  //parlist += GetFidCut()->GetFidCutParametersList()
2206 
2207  return new TObjString(parList) ;
2208 }
2209 
2210 //_________________________________________________________________________________
2212 //_________________________________________________________________________________
2213 void AliAnaCalorimeterQA::ExoticHistograms(Int_t absIdMax, Float_t ampMax,
2214  AliVCluster *clus, AliVCaloCells* cells)
2215 {
2216  if(ampMax < 0.01)
2217  {
2218  AliDebug(1,Form("Low amplitude energy %f",ampMax));
2219  return;
2220  }
2221 
2222  Float_t l0 = clus->GetM02();
2223  Float_t l1 = clus->GetM20();
2224  Float_t en = clus->E();
2225  Int_t nc = clus->GetNCells();
2226  Double_t tmax = clus->GetTOF()*1.e9-fConstantTimeShift; // recalibrated elsewhere
2227 
2228  Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
2229 
2230  if(en > 5)
2231  {
2232  fhExoL0ECross->Fill(eCrossFrac, l0, GetEventWeight());
2233  fhExoL1ECross->Fill(eCrossFrac, l1, GetEventWeight());
2234  }
2235 
2236  for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
2237  {
2238  for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
2239  {
2240  eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
2241 
2242  if(eCrossFrac > fExoECrossCuts[ie])
2243  {
2244  // Exotic
2245  fhExoL0 [ie][idt]->Fill(en,l0 , GetEventWeight());
2246  fhExoL1 [ie][idt]->Fill(en,l1 , GetEventWeight());
2247  fhExoTime [ie][idt]->Fill(en,tmax, GetEventWeight());
2248 
2249  if(en > 5)
2250  {
2251  fhExoL0NCell[ie][idt]->Fill(nc, l0, GetEventWeight());
2252  fhExoL1NCell[ie][idt]->Fill(nc, l1, GetEventWeight());
2253  }
2254 
2255  // Diff time, do for one cut in time
2256  if(idt == 0)
2257  {
2258  for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
2259  {
2260  Int_t absId = clus->GetCellsAbsId()[icell];
2261  Double_t time = cells->GetCellTime(absId);
2262  GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2263  time *= 1e9;
2264  time -= fConstantTimeShift;
2265 
2266  Float_t diff = tmax-time;
2267  fhExoDTime[ie]->Fill(en, diff, GetEventWeight());
2268  }
2269  }
2270  }
2271  else
2272  {
2273  fhExoECross[ie][idt]->Fill(en, eCrossFrac, GetEventWeight());
2274  fhExoNCell [ie][idt]->Fill(en, nc , GetEventWeight());
2275  }
2276  } // D time cut loop
2277  } // e cross cut loop
2278 }
2279 
2280 //___________________________________________________
2283 //___________________________________________________
2284 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
2285 {
2286  TList * outputContainer = new TList() ;
2287  outputContainer->SetName("QAHistos") ;
2288 
2289  // Init the number of modules, set in the class AliCalorimeterUtils
2290 
2291  fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
2292  if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
2293 
2294  // Histogram binning and ranges
2295 
2296  Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
2297  Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins(); Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax(); Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
2298  Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
2299  Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
2300  Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins(); Float_t massmax = GetHistogramRanges()->GetHistoMassMax(); Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
2301  Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins(); Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax(); Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
2302  Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
2303  Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
2304  Int_t ndRbins = GetHistogramRanges()->GetHistodRBins(); Float_t dRmax = GetHistogramRanges()->GetHistodRMax(); Float_t dRmin = GetHistogramRanges()->GetHistodRMin();
2305  Int_t ntimebins = GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
2306  Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins(); Int_t nclmax = GetHistogramRanges()->GetHistoNClustersMax(); Int_t nclmin = GetHistogramRanges()->GetHistoNClustersMin();
2307  Int_t ncebins = GetHistogramRanges()->GetHistoNCellsBins(); Int_t ncemax = GetHistogramRanges()->GetHistoNCellsMax(); Int_t ncemin = GetHistogramRanges()->GetHistoNCellsMin();
2308  Int_t nceclbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nceclmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nceclmin = GetHistogramRanges()->GetHistoNClusterCellMin();
2309  Int_t nvdistbins = GetHistogramRanges()->GetHistoVertexDistBins(); Float_t vdistmax = GetHistogramRanges()->GetHistoVertexDistMax(); Float_t vdistmin = GetHistogramRanges()->GetHistoVertexDistMin();
2310  Int_t rbins = GetHistogramRanges()->GetHistoRBins(); Float_t rmax = GetHistogramRanges()->GetHistoRMax(); Float_t rmin = GetHistogramRanges()->GetHistoRMin();
2311  Int_t xbins = GetHistogramRanges()->GetHistoXBins(); Float_t xmax = GetHistogramRanges()->GetHistoXMax(); Float_t xmin = GetHistogramRanges()->GetHistoXMin();
2312  Int_t ybins = GetHistogramRanges()->GetHistoYBins(); Float_t ymax = GetHistogramRanges()->GetHistoYMax(); Float_t ymin = GetHistogramRanges()->GetHistoYMin();
2313  Int_t zbins = GetHistogramRanges()->GetHistoZBins(); Float_t zmax = GetHistogramRanges()->GetHistoZMax(); Float_t zmin = GetHistogramRanges()->GetHistoZMin();
2314  Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
2315  Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
2316 
2317  Int_t nv0sbins = GetHistogramRanges()->GetHistoV0SignalBins(); Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax(); Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
2318  Int_t nv0mbins = GetHistogramRanges()->GetHistoV0MultiplicityBins(); Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax(); Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
2319  Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
2320 
2321  // Calorimeter grid
2322 
2323  // EMCAL
2324  fNMaxCols = 48;
2325  fNMaxRows = 24;
2326  fNRCU = 2 ;
2327  // PHOS
2328  if(GetCalorimeter()==kPHOS)
2329  {
2330  fNMaxCols = 56;
2331  fNMaxRows = 64;
2332  fNRCU = 4 ;
2333  }
2334 
2335 
2336  // Init histograms
2337 
2338  fhE = new TH1F ("hE","#it{E} reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
2339  fhE->SetXTitle("#it{E} (GeV)");
2340  outputContainer->Add(fhE);
2341 
2342  fhPt = new TH1F ("hPt","#it{p}_{T} reconstructed clusters", nptbins,ptmin,ptmax);
2343  fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2344  outputContainer->Add(fhPt);
2345 
2346  fhPhi = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
2347  fhPhi->SetXTitle("#phi (rad)");
2348  outputContainer->Add(fhPhi);
2349 
2350  fhEta = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
2351  fhEta->SetXTitle("#eta ");
2352  outputContainer->Add(fhEta);
2353 
2354  if(!fFillEBinAcceptanceHisto)
2355  {
2356  if(fFillAllTH3)
2357  {
2358  fhEtaPhiE = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
2359  netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2360  fhEtaPhiE->SetXTitle("#eta ");
2361  fhEtaPhiE->SetYTitle("#phi (rad)");
2362  fhEtaPhiE->SetZTitle("#it{E} (GeV) ");
2363  outputContainer->Add(fhEtaPhiE);
2364  }
2365  else
2366  {
2367  fhEtaPhi = new TH2F ("hEtaPhi","#eta vs #phi for #it{E} > 0.5 GeV, reconstructed clusters",
2368  netabins,etamin,etamax,nphibins,phimin,phimax);
2369  fhEtaPhi->SetXTitle("#eta ");
2370  fhEtaPhi->SetYTitle("#phi (rad)");
2371  outputContainer->Add(fhEtaPhi);
2372  }
2373  }
2374 
2375  fhClusterTimeEnergy = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
2376  nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2377  fhClusterTimeEnergy->SetXTitle("#it{E} (GeV) ");
2378  fhClusterTimeEnergy->SetYTitle("TOF (ns)");
2379  outputContainer->Add(fhClusterTimeEnergy);
2380 
2381  fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
2382  nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2383  fhClusterPairDiffTimeE->SetXTitle("#it{E}_{cluster} (GeV)");
2384  fhClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
2385  outputContainer->Add(fhClusterPairDiffTimeE);
2386 
2387  fhClusterPairDiffTimeESameMod = new TH2F("hClusterPairDiffTimeESameMod","cluster pair time difference vs E, only good clusters",
2388  nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2389  fhClusterPairDiffTimeESameMod->SetXTitle("#it{E}_{cluster} (GeV)");
2390  fhClusterPairDiffTimeESameMod->SetYTitle("#Delta #it{t} (ns)");
2391  outputContainer->Add(fhClusterPairDiffTimeESameMod);
2392 
2393  fhLambda0 = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
2394  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2395  fhLambda0->SetXTitle("#it{E}_{cluster}");
2396  fhLambda0->SetYTitle("#lambda^{2}_{0}");
2397  outputContainer->Add(fhLambda0);
2398 
2399  fhLambda1 = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E",
2400  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2401  fhLambda1->SetXTitle("#it{E}_{cluster}");
2402  fhLambda1->SetYTitle("#lambda^{2}_{1}");
2403  outputContainer->Add(fhLambda1);
2404 
2405  if(fStudyFECCorrelation)
2406  {
2407  fhLambda0MaxFECPair = new TH2F ("hLambda0MaxFECPair","#lambda^{2}_{0} vs E, FEC pair eta column",
2408  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2409  fhLambda0MaxFECPair->SetXTitle("#it{E}_{cluster}");
2410  fhLambda0MaxFECPair->SetYTitle("#lambda^{2}_{0}");
2411  outputContainer->Add(fhLambda0MaxFECPair);
2412 
2413  fhLambda1MaxFECPair = new TH2F ("hLambda1MaxFECPair","#lambda^{2}_{1} vs E, FEC pair eta column",
2414  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2415  fhLambda1MaxFECPair->SetXTitle("#it{E}_{cluster}");
2416  fhLambda1MaxFECPair->SetYTitle("#lambda^{2}_{1}");
2417  outputContainer->Add(fhLambda1MaxFECPair);
2418 
2419  fhLambda0Max2FECPair = new TH2F ("hLambda0Max2FECPair","#lambda^{2}_{0} vs E, FEC pair eta column, 2 highest cells corr.",
2420  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2421  fhLambda0Max2FECPair->SetXTitle("#it{E}_{cluster}");
2422  fhLambda0Max2FECPair->SetYTitle("#lambda^{2}_{0}");
2423  outputContainer->Add(fhLambda0Max2FECPair);
2424 
2425  fhLambda1Max2FECPair = new TH2F ("hLambda1Max2FECPair","#lambda^{2}_{1} vs E, FEC pair eta column, 2 highest cells corr.",
2426  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2427  fhLambda1Max2FECPair->SetXTitle("#it{E}_{cluster}");
2428  fhLambda1Max2FECPair->SetYTitle("#lambda^{2}_{1}");
2429  outputContainer->Add(fhLambda1Max2FECPair);
2430 
2431  fhLambda0MaxFECOdd = new TH2F ("hLambda0MaxFECOdd","#lambda^{2}_{0} vs E, FEC odd eta column",
2432  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2433  fhLambda0MaxFECOdd->SetXTitle("#it{E}_{cluster}");
2434  fhLambda0MaxFECOdd->SetYTitle("#lambda^{2}_{0}");
2435  outputContainer->Add(fhLambda0MaxFECOdd);
2436 
2437  fhLambda1MaxFECOdd = new TH2F ("hLambda1MaxFECOdd","#lambda^{2}_{1} vs E, FEC odd eta column",
2438  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2439  fhLambda1MaxFECOdd->SetXTitle("#it{E}_{cluster}");
2440  fhLambda1MaxFECOdd->SetYTitle("#lambda^{2}_{1}");
2441  outputContainer->Add(fhLambda1MaxFECOdd);
2442 
2443  fhLambda0Max2FECOdd = new TH2F ("hLambda0Max2FECOdd","#lambda^{2}_{0} vs E, FEC odd eta column, 2 highest cells corr.",
2444  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2445  fhLambda0Max2FECOdd->SetXTitle("#it{E}_{cluster}");
2446  fhLambda0Max2FECOdd->SetYTitle("#lambda^{2}_{0}");
2447  outputContainer->Add(fhLambda0Max2FECOdd);
2448 
2449  fhLambda1Max2FECOdd = new TH2F ("hLambda1Max2FECOdd","#lambda^{2}_{1} vs E, FEC odd eta column, 2 highest cells corr.",
2450  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2451  fhLambda1Max2FECOdd->SetXTitle("#it{E}_{cluster}");
2452  fhLambda1Max2FECOdd->SetYTitle("#lambda^{2}_{1}");
2453  outputContainer->Add(fhLambda1Max2FECOdd);
2454 
2455  fhLambda0MaxFECPairLargeNCells = new TH2F ("hLambda0MaxFECPair_LargeNCells","#lambda^{2}_{0} vs E, FEC pair eta column, ncell>5",
2456  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2457  fhLambda0MaxFECPairLargeNCells->SetXTitle("#it{E}_{cluster}");
2458  fhLambda0MaxFECPairLargeNCells->SetYTitle("#lambda^{2}_{0}");
2459  outputContainer->Add(fhLambda0MaxFECPairLargeNCells);
2460 
2461  fhLambda1MaxFECPairLargeNCells = new TH2F ("hLambda1MaxFECPair_LargeNCells","#lambda^{2}_{1} vs E, FEC pair eta column, ncell>5",
2462  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2463  fhLambda1MaxFECPairLargeNCells->SetXTitle("#it{E}_{cluster}");
2464  fhLambda1MaxFECPairLargeNCells->SetYTitle("#lambda^{2}_{1}");
2465  outputContainer->Add(fhLambda1MaxFECPairLargeNCells);
2466 
2467  fhLambda0Max2FECPairLargeNCells = new TH2F ("hLambda0Max2FECPair_LargeNCells","#lambda^{2}_{0} vs E, FEC pair eta column, 2 highest cells corr., ncell>5",
2468  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2469  fhLambda0Max2FECPairLargeNCells->SetXTitle("#it{E}_{cluster}");
2470  fhLambda0Max2FECPairLargeNCells->SetYTitle("#lambda^{2}_{0}");
2471  outputContainer->Add(fhLambda0Max2FECPairLargeNCells);
2472 
2473  fhLambda1Max2FECPairLargeNCells = new TH2F ("hLambda1Max2FECPair_LargeNCells","#lambda^{2}_{1} vs E, FEC pair eta column, 2 highest cells corr., ncell>5",
2474  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2475  fhLambda1Max2FECPairLargeNCells->SetXTitle("#it{E}_{cluster}");
2476  fhLambda1Max2FECPairLargeNCells->SetYTitle("#lambda^{2}_{1}");
2477  outputContainer->Add(fhLambda1Max2FECPairLargeNCells);
2478 
2479  fhLambda0MaxFECOddLargeNCells = new TH2F ("hLambda0MaxFECOdd_LargeNCells","#lambda^{2}_{0} vs E, FEC odd eta column, ncell>5",
2480  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2481  fhLambda0MaxFECOddLargeNCells->SetXTitle("#it{E}_{cluster}");
2482  fhLambda0MaxFECOddLargeNCells->SetYTitle("#lambda^{2}_{0}");
2483  outputContainer->Add(fhLambda0MaxFECOddLargeNCells);
2484 
2485  fhLambda1MaxFECOddLargeNCells = new TH2F ("hLambda1MaxFECOdd_LargeNCells","#lambda^{2}_{1} vs E, FEC odd eta column, ncell>5",
2486  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2487  fhLambda1MaxFECOddLargeNCells->SetXTitle("#it{E}_{cluster}");
2488  fhLambda1MaxFECOddLargeNCells->SetYTitle("#lambda^{2}_{1}");
2489  outputContainer->Add(fhLambda1MaxFECOddLargeNCells);
2490 
2491  fhLambda0Max2FECOddLargeNCells = new TH2F ("hLambda0Max2FECOdd_LargeNCells","#lambda^{2}_{0} vs E, FEC odd eta column, 2 highest cells corr., ncell>5",
2492  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2493  fhLambda0Max2FECOddLargeNCells->SetXTitle("#it{E}_{cluster}");
2494  fhLambda0Max2FECOddLargeNCells->SetYTitle("#lambda^{2}_{0}");
2495  outputContainer->Add(fhLambda0Max2FECOddLargeNCells);
2496 
2497  fhLambda1Max2FECOddLargeNCells = new TH2F ("hLambda1Max2FECOdd_LargeNCells","#lambda^{2}_{1} vs E, FEC odd eta column, 2 highest cells corr., ncell>5",
2498  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2499  fhLambda1Max2FECOddLargeNCells->SetXTitle("#it{E}_{cluster}");
2500  fhLambda1Max2FECOddLargeNCells->SetYTitle("#lambda^{2}_{1}");
2501  outputContainer->Add(fhLambda1Max2FECOddLargeNCells);
2502 
2503 
2504  TString titleMaxFEC[] =
2505  {"No correlation with main cell", "Main cell corr. with near cell",
2506  "Main cell corr. with far cell", "Main cell corr. with near and far"};
2507 
2508  TString titleFEC[] =
2509  {"No main cell corr., no other inside cluster" , "No main cell corr., 1 cell correl in cluster" , "No main cell corr., >=2 cell correl in cluster",
2510  "Main cell corr. & no other inside cluster" , "Main cell corr. & 1 cell corr. inside cluster", "Main cell corr. & >=2 cell corr. inside cluster"};
2511 
2512  fhEtaPhiFECCorrControl = new TH2F ("hEtaPhiFECCorrControl","#eta vs #phi for FEC correlation selected clusters",
2513  netabins,etamin,etamax,nphibins,phimin,phimax);
2514  fhEtaPhiFECCorrControl->SetXTitle("#eta ");
2515  fhEtaPhiFECCorrControl->SetYTitle("#phi (rad)");
2516  outputContainer->Add(fhEtaPhiFECCorrControl);
2517 
2518  for(Int_t i = 0; i < 4; i++)
2519  {
2520  fhLambda0MaxFECCorrel[i] = new TH2F
2521  (Form("hLambda0MaxFECCorrel_Case%d",i),
2522  Form("#lambda^{2}_{0} vs E, Max FEC correl %s",titleMaxFEC[i].Data()),
2523  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2524  fhLambda0MaxFECCorrel[i]->SetXTitle("#it{E}_{cluster}");
2525  fhLambda0MaxFECCorrel[i]->SetYTitle("#lambda^{2}_{0}");
2526  outputContainer->Add(fhLambda0MaxFECCorrel[i]);
2527 
2528  fhLambda1MaxFECCorrel[i] = new TH2F
2529  (Form("hLambda1MaxFECCorrel_Case%d",i),
2530  Form("#lambda^{2}_{1} vs E for bad cluster, Max FEC correl %s",titleMaxFEC[i].Data()),
2531  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2532  fhLambda1MaxFECCorrel[i]->SetXTitle("#it{E}_{cluster}");
2533  fhLambda1MaxFECCorrel[i]->SetYTitle("#lambda^{2}_{1}");
2534  outputContainer->Add(fhLambda1MaxFECCorrel[i]);
2535 
2536  fhLambda0MaxFECCorrelLargeNCells[i] = new TH2F
2537  (Form("hLambda0MaxFECCorrel_Case%d_LargeNCells",i),
2538  Form("#lambda^{2}_{0} vs E, Max FEC correl %s, ncells>5",titleMaxFEC[i].Data()),
2539  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2540  fhLambda0MaxFECCorrelLargeNCells[i]->SetXTitle("#it{E}_{cluster}");
2541  fhLambda0MaxFECCorrelLargeNCells[i]->SetYTitle("#lambda^{2}_{0}");
2542  outputContainer->Add(fhLambda0MaxFECCorrelLargeNCells[i]);
2543 
2544  fhLambda1MaxFECCorrelLargeNCells[i] = new TH2F
2545  (Form("hLambda1MaxFECCorrel_Case%d_LargeNCells",i),
2546  Form("#lambda^{2}_{1} vs E for bad cluster, Max FEC correl %s, ncells>5",titleMaxFEC[i].Data()),
2547  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2548  fhLambda1MaxFECCorrelLargeNCells[i]->SetXTitle("#it{E}_{cluster}");
2549  fhLambda1MaxFECCorrelLargeNCells[i]->SetYTitle("#lambda^{2}_{1}");
2550  outputContainer->Add(fhLambda1MaxFECCorrelLargeNCells[i]);
2551  }
2552 
2553  for(Int_t i = 0; i < 6; i++)
2554  {
2555  fhLambda0FECCorrel[i] = new TH2F
2556  (Form("hLambda0FECCorrel_Case%d",i),
2557  Form("#lambda^{2}_{0} vs E, %s",titleFEC[i].Data()),
2558  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2559  fhLambda0FECCorrel[i]->SetXTitle("#it{E}_{cluster}");
2560  fhLambda0FECCorrel[i]->SetYTitle("#lambda^{2}_{0}");
2561  outputContainer->Add(fhLambda0FECCorrel[i]);
2562 
2563  fhLambda1FECCorrel[i] = new TH2F
2564  (Form("hLambda1FECCorrel_Case%d",i),
2565  Form("#lambda^{2}_{1} vs E for bad cluster, %s",titleFEC[i].Data()),
2566  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2567  fhLambda1FECCorrel[i]->SetXTitle("#it{E}_{cluster}");
2568  fhLambda1FECCorrel[i]->SetYTitle("#lambda^{2}_{1}");
2569  outputContainer->Add(fhLambda1FECCorrel[i]);
2570 
2571  fhLambda0FECCorrelLargeNCells[i] = new TH2F
2572  (Form("hLambda0FECCorrel_Case%d_LargeNCells",i),
2573  Form("#lambda^{2}_{0} vs E, %s, ncells>5",titleFEC[i].Data()),
2574  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2575  fhLambda0FECCorrelLargeNCells[i]->SetXTitle("#it{E}_{cluster}");
2576  fhLambda0FECCorrelLargeNCells[i]->SetYTitle("#lambda^{2}_{0}");
2577  outputContainer->Add(fhLambda0FECCorrelLargeNCells[i]);
2578 
2579  fhLambda1FECCorrelLargeNCells[i] = new TH2F
2580  (Form("hLambda1FECCorrel_Case%d_LargeNCells",i),
2581  Form("#lambda^{2}_{1} vs E for bad cluster, %s, ncells>5",titleFEC[i].Data()),
2582  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2583  fhLambda1FECCorrelLargeNCells[i]->SetXTitle("#it{E}_{cluster}");
2584  fhLambda1FECCorrelLargeNCells[i]->SetYTitle("#lambda^{2}_{1}");
2585  outputContainer->Add(fhLambda1FECCorrelLargeNCells[i]);
2586  }
2587 
2588  }
2589 // fhDispersion = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
2590 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2591 // fhDispersion->SetXTitle("#it{E}_{cluster}");
2592 // fhDispersion->SetYTitle("Dispersion");
2593 // outputContainer->Add(fhDispersion);
2594 
2595  fhClusterMaxCellCloseCellRatio = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
2596  nptbins,ptmin,ptmax, 100,0,1.);
2597  fhClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2598  fhClusterMaxCellCloseCellRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cell max}");
2599  outputContainer->Add(fhClusterMaxCellCloseCellRatio);
2600 
2601  fhClusterMaxCellCloseCellDiff = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
2602  nptbins,ptmin,ptmax, 500,0,100.);
2603  fhClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
2604  fhClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max}-#it{E}_{cell i} (GeV)");
2605  outputContainer->Add(fhClusterMaxCellCloseCellDiff);
2606 
2607  fhClusterMaxCellDiff = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
2608  nptbins,ptmin,ptmax, 500,0,1.);
2609  fhClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
2610  fhClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
2611  outputContainer->Add(fhClusterMaxCellDiff);
2612 
2613  if(fStudyBadClusters)
2614  {
2615  fhClusterMaxCellDiffNoCut = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
2616  nptbins,ptmin,ptmax, 500,0,1.);
2617  fhClusterMaxCellDiffNoCut->SetXTitle("#it{E}_{cluster} (GeV) ");
2618  fhClusterMaxCellDiffNoCut->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
2619  outputContainer->Add(fhClusterMaxCellDiffNoCut);
2620  }
2621 
2622  fhClusterMaxCellECross = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
2623  nptbins,ptmin,ptmax, 400,-1,1.);
2624  fhClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
2625  fhClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
2626  outputContainer->Add(fhClusterMaxCellECross);
2627 
2628  if(fStudyBadClusters)
2629  {
2630  fhNCellsPerClusterNoCut = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
2631  nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2632  fhNCellsPerClusterNoCut->SetXTitle("#it{E} (GeV)");
2633  fhNCellsPerClusterNoCut->SetYTitle("#it{n}_{cells}");
2634  outputContainer->Add(fhNCellsPerClusterNoCut);
2635  }
2636 
2637  fhNCellsPerCluster = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2638  fhNCellsPerCluster->SetXTitle("#it{E} (GeV)");
2639  fhNCellsPerCluster->SetYTitle("#it{n}_{cells}");
2640  outputContainer->Add(fhNCellsPerCluster);
2641 
2642  if(fStudyFECCorrelation)
2643  {
2644  fhNCellsPerClusterWithWeight = new TH2F
2645  ("hNCellsPerClusterWithWeight","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2646  fhNCellsPerClusterWithWeight->SetXTitle("#it{E} (GeV)");
2647  fhNCellsPerClusterWithWeight->SetYTitle("#it{n}_{cells}");
2648  outputContainer->Add(fhNCellsPerClusterWithWeight);
2649 
2650  fhNCellsPerClusterRatioWithWeight = new TH2F
2651  ("hNCellsPerClusterRatioWithWeight","# cells per cluster vs energy",nptbins,ptmin,ptmax, 100,0,1);
2652  fhNCellsPerClusterRatioWithWeight->SetXTitle("#it{E} (GeV)");
2653  fhNCellsPerClusterRatioWithWeight->SetYTitle("#it{n}_{cells}");
2654  outputContainer->Add(fhNCellsPerClusterRatioWithWeight);
2655  }
2656 
2657  if(fStudyBadClusters)
2658  {
2659  fhNCellsPerClusterWeirdModNoCut = new TH2F ("hNCellsPerClusterWeirdNoCutMod","# cells per cluster vs energy, E > 100, no bad clusters cut, per SM",
2660  nceclbins,nceclmin,nceclmax,fNModules,0,fNModules);
2661  fhNCellsPerClusterWeirdModNoCut->SetYTitle("SM number");
2662  fhNCellsPerClusterWeirdModNoCut->SetXTitle("#it{n}_{cells}");
2663  outputContainer->Add(fhNCellsPerClusterWeirdModNoCut);
2664  }
2665 
2666  fhNCellsPerClusterWeirdMod = new TH2F ("hNCellsPerClusterWeirdMod","# cells per cluster, E > 100 GeV, per SM",
2667  nceclbins*2,nceclmin,nceclmax*2,fNModules,0,fNModules);
2668  fhNCellsPerClusterWeirdMod->SetYTitle("SM number");
2669  fhNCellsPerClusterWeirdMod->SetXTitle("#it{n}_{cells}");
2670  outputContainer->Add(fhNCellsPerClusterWeirdMod);
2671 
2672  fhNClusters = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
2673  fhNClusters->SetXTitle("#it{n}_{clusters}");
2674  outputContainer->Add(fhNClusters);
2675 
2676  if(fStudyBadClusters)
2677  {
2678  fhBadClusterEnergy = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
2679  fhBadClusterEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
2680  outputContainer->Add(fhBadClusterEnergy);
2681 
2682  fhBadClusterEtaPhi = new TH2F ("hBadClusterEtaPhi","Bad cluster, #eta vs #phi, #it{E} > 0.5 GeV",
2683  netabins,etamin,etamax,nphibins,phimin,phimax);
2684  fhBadClusterEtaPhi->SetXTitle("#eta ");
2685  fhBadClusterEtaPhi->SetXTitle("#phi (rad) ");
2686  outputContainer->Add(fhBadClusterEtaPhi);
2687 
2688  fhBadClusterLambda0 = new TH2F ("hBadClusterLambda0","Bad cluster,shower shape, #lambda^{2}_{0} vs E",
2689  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2690  fhBadClusterLambda0->SetXTitle("#it{E}_{cluster}");
2691  fhBadClusterLambda0->SetYTitle("#lambda^{2}_{0}");
2692  outputContainer->Add(fhBadClusterLambda0);
2693 
2694  fhBadClusterLambda1 = new TH2F ("hBadClusterLambda1","Bad cluster,shower shape, #lambda^{2}_{1} vs E for bad cluster ",
2695  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2696  fhBadClusterLambda1->SetXTitle("#it{E}_{cluster}");
2697  fhBadClusterLambda1->SetYTitle("#lambda^{2}_{1}");
2698  outputContainer->Add(fhBadClusterLambda1);
2699 
2700  fhBadClusterMaxCellCloseCellRatio = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
2701  nptbins,ptmin,ptmax, 100,0,1.);
2702  fhBadClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2703  fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
2704  outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
2705 
2706  fhBadClusterMaxCellCloseCellDiff = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
2707  nptbins,ptmin,ptmax, 500,0,100);
2708  fhBadClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
2709  fhBadClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max} - #it{E}_{cell i} (GeV)");
2710  outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
2711 
2712  fhBadClusterMaxCellDiff = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
2713  nptbins,ptmin,ptmax, 500,0,1.);
2714  fhBadClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
2715  fhBadClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max}) / #it{E}_{cluster}");
2716  outputContainer->Add(fhBadClusterMaxCellDiff);
2717 
2718  fhBadClusterTimeEnergy = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
2719  nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2720  fhBadClusterTimeEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
2721  fhBadClusterTimeEnergy->SetYTitle("#it{t} (ns)");
2722  outputContainer->Add(fhBadClusterTimeEnergy);
2723 
2724  fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2725  fhBadClusterPairDiffTimeE->SetXTitle("#it{E}_{bad cluster} (GeV)");
2726  fhBadClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
2727  outputContainer->Add(fhBadClusterPairDiffTimeE);
2728 
2729  fhBadClusterMaxCellECross = new TH2F ("hBadClusterMaxCellECross","1 - #it{E}_{+} around max energy cell / max energy cell vs cluster energy, bad clusters",
2730  nptbins,ptmin,ptmax, 400,-1,1.);
2731  fhBadClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
2732  fhBadClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
2733  outputContainer->Add(fhBadClusterMaxCellECross);
2734 
2735  if(fFillAllCellTimeHisto)
2736  {
2737  fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","#it{t}_{cell max}-#it{t}_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2738  fhBadCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
2739  fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max - i} (ns)");
2740  outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
2741 
2742  fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","#it{t}_{cell max}-#it{t}_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2743  fhBadClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
2744  fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
2745  outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
2746 
2747  fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","#it{t}_{cell max}-#it{t}_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2748  fhBadClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
2749  fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
2750  outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
2751 
2752  }
2753 
2754  }
2755 
2756  if(fStudyExotic)
2757  {
2758  fhExoL0ECross = new TH2F("hExoL0_ECross",
2759  "#lambda^{2}_{0} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
2760  400,0,1,ssbins,ssmin,ssmax);
2761  fhExoL0ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
2762  fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
2763  outputContainer->Add(fhExoL0ECross) ;
2764 
2765  fhExoL1ECross = new TH2F("hExoL1_ECross",
2766  "#lambda^{2}_{1} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
2767  400,0,1,ssbins,ssmin,ssmax);
2768  fhExoL1ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
2769  fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
2770  outputContainer->Add(fhExoL1ECross) ;
2771 
2772  for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
2773  {
2774 
2775  fhExoDTime[ie] = new TH2F(Form("hExoDTime_ECross%d",ie),
2776  Form("#Delta time = t_{max}-t_{cells} vs #it{E}_{cluster} for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f",fExoECrossCuts[ie]),
2777  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
2778  fhExoDTime[ie] ->SetYTitle("#Delta #it{t} (ns)");
2779  fhExoDTime[ie] ->SetXTitle("#it{E} (GeV)");
2780  outputContainer->Add(fhExoDTime[ie]) ;
2781 
2782  for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
2783  {
2784  fhExoNCell[ie][idt] = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
2785  Form("N cells per cluster vs E cluster, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2786  nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
2787  fhExoNCell[ie][idt] ->SetYTitle("#it{n}_cells");
2788  fhExoNCell[ie][idt] ->SetXTitle("#it{E} (GeV)");
2789  outputContainer->Add(fhExoNCell[ie][idt]) ;
2790 
2791  fhExoL0 [ie][idt] = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
2792  Form("#lambda^{2}_{0} vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2793  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2794  fhExoL0 [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
2795  fhExoL0 [ie][idt] ->SetXTitle("#it{E} (GeV)");
2796  outputContainer->Add(fhExoL0[ie][idt]) ;
2797 
2798  fhExoL1 [ie][idt] = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
2799  Form("#lambda^{2}_{1} vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2800  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2801  fhExoL1 [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
2802  fhExoL1 [ie][idt] ->SetXTitle("#it{E} (GeV)");
2803  outputContainer->Add(fhExoL1[ie][idt]) ;
2804 
2805  fhExoECross[ie][idt] = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
2806  Form("#it{E} cross for cells vs E cell, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2807  nptbins,ptmin,ptmax,400,0,1);
2808  fhExoECross[ie][idt] ->SetYTitle("1-#it{E}_{+}/#it{E}_{cell max}");
2809  fhExoECross[ie][idt] ->SetXTitle("#it{E}_{cell} (GeV)");
2810  outputContainer->Add(fhExoECross[ie][idt]) ;
2811 
2812  fhExoTime [ie][idt] = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
2813  Form("Time of cluster (max cell) vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2814  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2815  fhExoTime [ie][idt] ->SetYTitle("#it{t}_{max} (ns)");
2816  fhExoTime [ie][idt] ->SetXTitle("#it{E} (GeV)");
2817  outputContainer->Add(fhExoTime[ie][idt]) ;
2818 
2819  fhExoL0NCell[ie][idt] = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
2820  Form("#lambda^{2}_{0} vs N cells per clusters for E > 5 GeV, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2821  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2822  fhExoL0NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
2823  fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
2824  outputContainer->Add(fhExoL0NCell[ie][idt]) ;
2825 
2826  fhExoL1NCell[ie][idt] = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
2827  Form("#lambda^{2}_{1} vs N cells per clusters for E > 5 GeV, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
2828  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2829  fhExoL1NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
2830  fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
2831  outputContainer->Add(fhExoL1NCell[ie][idt]) ;
2832  }
2833  }
2834  }
2835 
2836  // Cluster size in terms of cells
2837 
2838  if(fStudyClustersAsymmetry)
2839  {
2840  fhDeltaIEtaDeltaIPhiE0[0] = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2841  50,0,50,50,0,50);
2842  fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
2843  fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
2844  outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
2845 
2846  fhDeltaIEtaDeltaIPhiE2[0] = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2847  50,0,50,50,0,50);
2848  fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
2849  fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
2850  outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
2851 
2852  fhDeltaIEtaDeltaIPhiE6[0] = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2853  50,0,50,50,0,50);
2854  fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
2855  fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
2856  outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
2857 
2858  fhDeltaIA[0] = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
2859  nptbins,ptmin,ptmax,21,-1.05,1.05);
2860  fhDeltaIA[0]->SetXTitle("#it{E}_{cluster}");
2861  fhDeltaIA[0]->SetYTitle("#it{A}_{cell in cluster}");
2862  outputContainer->Add(fhDeltaIA[0]);
2863 
2864  fhDeltaIAL0[0] = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
2865  ssbins,ssmin,ssmax,21,-1.05,1.05);
2866  fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
2867  fhDeltaIAL0[0]->SetYTitle("#it{A}_{cell in cluster}");
2868  outputContainer->Add(fhDeltaIAL0[0]);
2869 
2870  fhDeltaIAL1[0] = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
2871  ssbins,ssmin,ssmax,21,-1.05,1.05);
2872  fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
2873  fhDeltaIAL1[0]->SetYTitle("#it{A}_{cell in cluster}");
2874  outputContainer->Add(fhDeltaIAL1[0]);
2875 
2876  fhDeltaIANCells[0] = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
2877  nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
2878  fhDeltaIANCells[0]->SetXTitle("#it{n}_{cell in cluster}");
2879  fhDeltaIANCells[0]->SetYTitle("#it{A}_{cell in cluster}");
2880  outputContainer->Add(fhDeltaIANCells[0]);
2881 
2882 
2883  fhDeltaIEtaDeltaIPhiE0[1] = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3, matched with track",
2884  50,0,50,50,0,50);
2885  fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
2886  fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
2887  outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
2888 
2889  fhDeltaIEtaDeltaIPhiE2[1] = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3, matched with track",
2890  50,0,50,50,0,50);
2891  fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
2892  fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
2893  outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
2894 
2895  fhDeltaIEtaDeltaIPhiE6[1] = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3, matched with track",
2896  50,0,50,50,0,50);
2897  fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
2898  fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
2899  outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
2900 
2901  fhDeltaIA[1] = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
2902  nptbins,ptmin,ptmax,21,-1.05,1.05);
2903  fhDeltaIA[1]->SetXTitle("#it{E}_{cluster}");
2904  fhDeltaIA[1]->SetYTitle("#it{A}_{cell in cluster}");
2905  outputContainer->Add(fhDeltaIA[1]);
2906 
2907  fhDeltaIAL0[1] = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
2908  ssbins,ssmin,ssmax,21,-1.05,1.05);
2909  fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
2910  fhDeltaIAL0[1]->SetYTitle("#it{A}_{cell in cluster}");
2911  outputContainer->Add(fhDeltaIAL0[1]);
2912 
2913  fhDeltaIAL1[1] = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
2914  ssbins,ssmin,ssmax,21,-1.05,1.05);
2915  fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
2916  fhDeltaIAL1[1]->SetYTitle("#it{A}_{cell in cluster}");
2917  outputContainer->Add(fhDeltaIAL1[1]);
2918 
2919  fhDeltaIANCells[1] = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
2920  nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
2921  fhDeltaIANCells[1]->SetXTitle("#it{n}_{cell in cluster}");
2922  fhDeltaIANCells[1]->SetYTitle("#it{A}_{cell in cluster}");
2923  outputContainer->Add(fhDeltaIANCells[1]);
2924 
2925  if(IsDataMC()){
2926  TString particle[]={"Photon","Electron","Conversion","Hadron"};
2927  for (Int_t iPart = 0; iPart < 4; iPart++) {
2928 
2929  fhDeltaIAMC[iPart] = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
2930  nptbins,ptmin,ptmax,21,-1.05,1.05);
2931  fhDeltaIAMC[iPart]->SetXTitle("#it{E}_{cluster}");
2932  fhDeltaIAMC[iPart]->SetYTitle("#it{A}_{cell in cluster}");
2933  outputContainer->Add(fhDeltaIAMC[iPart]);
2934  }
2935  }
2936 
2937  if(fStudyBadClusters)
2938  {
2939  fhBadClusterDeltaIEtaDeltaIPhiE0 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2940  50,0,50,50,0,50);
2941  fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
2942  fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
2943  outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
2944 
2945  fhBadClusterDeltaIEtaDeltaIPhiE2 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2946  50,0,50,50,0,50);
2947  fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
2948  fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
2949  outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
2950 
2951  fhBadClusterDeltaIEtaDeltaIPhiE6 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2952  50,0,50,50,0,50);
2953  fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
2954  fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
2955  outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
2956 
2957  fhBadClusterDeltaIA = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
2958  nptbins,ptmin,ptmax,21,-1.05,1.05);
2959  fhBadClusterDeltaIA->SetXTitle("#it{E}_{cluster}");
2960  fhBadClusterDeltaIA->SetYTitle("#it{A}_{cell in cluster}");
2961  outputContainer->Add(fhBadClusterDeltaIA);
2962  }
2963  }
2964 
2965  if(fStudyWeight)
2966  {
2967  fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2968  nptbins,ptmin,ptmax, 100,0,1.);
2969  fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2970  fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
2971  outputContainer->Add(fhECellClusterRatio);
2972 
2973  fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
2974  nptbins,ptmin,ptmax, 100,-10,0);
2975  fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2976  fhECellClusterLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{cluster})");
2977  outputContainer->Add(fhECellClusterLogRatio);
2978 
2979  fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2980  nptbins,ptmin,ptmax, 100,0,1.);
2981  fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2982  fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
2983  outputContainer->Add(fhEMaxCellClusterRatio);
2984 
2985  fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
2986  nptbins,ptmin,ptmax, 100,-10,0);
2987  fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2988  fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
2989  outputContainer->Add(fhEMaxCellClusterLogRatio);
2990 
2991  fhECellTotalRatio = new TH2F ("hECellTotalRatio"," cell energy / sum all energy vs all energy",
2992  nptbins*2,ptmin,ptmax*2, 100,0,1.);
2993  fhECellTotalRatio->SetXTitle("#it{E}_{total} (GeV) ");
2994  fhECellTotalRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{total}");
2995  outputContainer->Add(fhECellTotalRatio);
2996 
2997  fhECellTotalLogRatio = new TH2F ("hECellTotalLogRatio"," Log(cell energy / sum all energy) vs all energy",
2998  nptbins*2,ptmin,ptmax*2, 100,-10,0);
2999  fhECellTotalLogRatio->SetXTitle("#it{E}_{total} (GeV) ");
3000  fhECellTotalLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{total})");
3001  outputContainer->Add(fhECellTotalLogRatio);
3002 
3003  fhECellTotalRatioMod = new TH2F*[fNModules];
3004  fhECellTotalLogRatioMod = new TH2F*[fNModules];
3005 
3006  for(Int_t imod = 0; imod < fNModules; imod++)
3007  {
3008  fhECellTotalRatioMod[imod] = new TH2F (Form("hECellTotalRatio_Mod%d",imod),
3009  Form("#cell energy / sum all energy vs all energy in Module %d",imod),
3010  nptbins*2,ptmin,ptmax*2, 100,0,1.);
3011  fhECellTotalRatioMod[imod]->SetXTitle("#it{E} (GeV)");
3012  fhECellTotalRatioMod[imod]->SetYTitle("#it{n}_{cells}");
3013  outputContainer->Add(fhECellTotalRatioMod[imod]);
3014 
3015  fhECellTotalLogRatioMod[imod] = new TH2F (Form("hECellTotalLogRatio_Mod%d",imod),
3016  Form("Log(cell energy / sum all energy) vs all energy in Module %d",imod),
3017  nptbins*2,ptmin,ptmax*2, 100,-10,0);
3018  fhECellTotalLogRatioMod[imod]->SetXTitle("#it{E} (GeV)");
3019  fhECellTotalLogRatioMod[imod]->SetYTitle("#it{n}_{cells}");
3020  outputContainer->Add(fhECellTotalLogRatioMod[imod]);
3021  }
3022 
3023  // To be done properly with setters and data members ...
3024  Float_t cellEmin [] = {0.05,0.1,0.15,0.2};
3025  Float_t cellTmin [] = {50.,100.,1000000.};
3026 
3027  for(Int_t iw = 0; iw < 12; iw++)
3028  {
3029  Float_t w0 = 4+0.05*iw; // 3+0.25*iw;
3030  for(Int_t iEmin = 0; iEmin < 4; iEmin++)
3031  {
3032  for(Int_t iTmin = 0; iTmin < 3; iTmin++)
3033  {
3034  fhLambda0ForW0AndCellCuts[iw][iEmin][iTmin] = new TH2F (Form("hLambda0ForW0%d_CellEMin%d_TimeMax%d",iw,iEmin,iTmin),
3035  Form("#lambda^{2}_{0} vs E, w0=%1.2f, cell E>%2.2f MeV, |t|<%2.0f ns",
3036  w0, cellEmin[iEmin], cellTmin[iTmin]),
3037  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3038  fhLambda0ForW0AndCellCuts[iw][iEmin][iTmin]->SetXTitle("#it{E}_{cluster}");
3039  fhLambda0ForW0AndCellCuts[iw][iEmin][iTmin]->SetYTitle("#lambda^{2}_{0}");
3040  outputContainer->Add(fhLambda0ForW0AndCellCuts[iw][iEmin][iTmin]);
3041 
3042 
3043 // fhLambda1ForW0AndCellCuts[iw][iEmin][iTmin] = new TH2F (Form("hLambda1ForW0%d_CellEMin%d_TimeMax%d",iw,iEmin,iTmin),
3044 // Form("#lambda^{2}_{1} vs E, w0=%1.2f, cell E>%2.2f MeV, |t|<%2.0f ns"",
3045 // w0, cellEmin[iEmin], cellTmin[iTmin]),
3046 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3047 // fhLambda1ForW0AndCellCuts[iw][iEmin][iTmin]->SetXTitle("#it{E}_{cluster}");
3048 // fhLambda1ForW0AndCellCuts[iw][iEmin][iTmin]->SetYTitle("#lambda^{2}_{1}");
3049 // outputContainer->Add(fhLambda1ForW0AndCellCuts[iw][iEmin][iTmin]);
3050 
3051  fhLambda0ForW0AndCellCutsEta0[iw][iEmin][iTmin] = new TH2F (Form("hLambda0ForW0%d_CellEMin%d_TimeMax%d_Eta0",iw,iEmin,iTmin),
3052  Form("#lambda^{2}_{0} vs E, w0=%1.2f, cell E>%2.2f MeV, |t|<%2.0f ns, |#eta| < 0.15",
3053  w0, cellEmin[iEmin], cellTmin[iTmin]),
3054  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3055  fhLambda0ForW0AndCellCutsEta0[iw][iEmin][iTmin]->SetXTitle("#it{E}_{cluster}");
3056  fhLambda0ForW0AndCellCutsEta0[iw][iEmin][iTmin]->SetYTitle("#lambda^{2}_{0}");
3057  outputContainer->Add(fhLambda0ForW0AndCellCutsEta0[iw][iEmin][iTmin]);
3058  }
3059  }
3060 
3061  if(IsDataMC())
3062  {
3063  TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
3064  for(Int_t imc = 0; imc < 5; imc++)
3065  {
3066  fhLambda0ForW0MC[iw][imc] = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
3067  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
3068  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3069  fhLambda0ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
3070  fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
3071  outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
3072 
3073 // fhLambda1ForW0MC[iw][imc] = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
3074 // Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
3075 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3076 // fhLambda1ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
3077 // fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
3078 // outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
3079  }
3080  }
3081  }
3082  }
3083 
3084  // Track Matching
3085  if(fFillAllTMHisto)
3086  {
3087  Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
3088  Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
3089  Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
3090  Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
3091  Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
3092  Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
3093 
3094  fhTrackMatchedDEtaNeg = new TH2F("hTrackMatchedDEtaNeg","d#eta of cluster-negative track vs cluster energy",
3095  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
3096  fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
3097  fhTrackMatchedDEtaNeg->SetXTitle("#it{E}_{cluster} (GeV)");
3098 
3099  fhTrackMatchedDPhiNeg = new TH2F("hTrackMatchedDPhiNeg","d#phi of cluster-negative track vs cluster energy",
3100  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
3101  fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
3102  fhTrackMatchedDPhiNeg->SetXTitle("#it{E}_{cluster} (GeV)");
3103 
3104  fhTrackMatchedDEtaDPhiNeg = new TH2F("hTrackMatchedDEtaDPhiNeg","d#eta vs d#phi of cluster- negative track vs cluster energy",
3105  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
3106  fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
3107  fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
3108 
3109  fhTrackMatchedDEtaPos = new TH2F("hTrackMatchedDEtaPos","d#eta of cluster-positive track vs cluster energy",
3110  nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
3111  fhTrackMatchedDEtaPos->SetYTitle("d#eta");
3112  fhTrackMatchedDEtaPos->SetXTitle("#it{E}_{cluster} (GeV)");
3113 
3114  fhTrackMatchedDPhiPos = new TH2F("hTrackMatchedDPhiPos","d#phi of cluster-positive track vs cluster energy",
3115  nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
3116  fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
3117  fhTrackMatchedDPhiPos->SetXTitle("#it{E}_{cluster} (GeV)");
3118 
3119  fhTrackMatchedDEtaDPhiPos = new TH2F("hTrackMatchedDEtaDPhiPos","d#eta vs d#phi of cluster-positive track vs cluster energy",
3120  nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
3121  fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
3122  fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
3123 
3124 
3125  fhTrackMatchedDEtaNegMod = new TH2F("hTrackMatchedDEtaNegPerModule","d#eta of cluster-negative track vs module, E > 0.5 GeV",
3126  nresetabins,resetamin,resetamax,fNModules,0,fNModules);
3127  fhTrackMatchedDEtaNegMod->SetXTitle("d#eta");
3128  fhTrackMatchedDEtaNegMod->SetYTitle("Module");
3129 
3130  fhTrackMatchedDPhiNegMod = new TH2F("hTrackMatchedDPhiNegPerModule","d#phi of cluster-negative track vs module, E > 0.5 GeV",
3131  nresetabins,resetamin,resetamax,fNModules,0,fNModules);
3132  fhTrackMatchedDPhiNegMod->SetXTitle("d#phi (rad)");
3133  fhTrackMatchedDPhiNegMod->SetYTitle("Module");
3134 
3135  fhTrackMatchedDEtaPosMod = new TH2F("hTrackMatchedDEtaPosPerModule","d#eta of cluster-positive track vs module, E > 0.5 GeV",
3136  nresetabins,resetamin,resetamax,fNModules,0,fNModules);
3137  fhTrackMatchedDEtaPosMod->SetXTitle("d#eta");
3138  fhTrackMatchedDEtaPosMod->SetYTitle("Module");
3139 
3140  fhTrackMatchedDPhiPosMod = new TH2F("hTrackMatchedDPhiPosPerModule","d#phi of cluster-positive track vs module, E > 0.5 GeV",
3141  nresetabins,resetamin,resetamax,fNModules,0,fNModules);
3142  fhTrackMatchedDPhiPosMod->SetXTitle("d#phi (rad)");
3143  fhTrackMatchedDPhiPosMod->SetYTitle("Module");
3144 
3145  outputContainer->Add(fhTrackMatchedDEtaNeg) ;
3146  outputContainer->Add(fhTrackMatchedDPhiNeg) ;
3147  outputContainer->Add(fhTrackMatchedDEtaPos) ;
3148  outputContainer->Add(fhTrackMatchedDPhiPos) ;
3149  outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
3150  outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
3151 
3152  outputContainer->Add(fhTrackMatchedDEtaNegMod) ;
3153  outputContainer->Add(fhTrackMatchedDPhiNegMod) ;
3154  outputContainer->Add(fhTrackMatchedDEtaPosMod) ;
3155  outputContainer->Add(fhTrackMatchedDPhiPosMod) ;
3156 
3157  fhECharged = new TH1F ("hECharged","#it{E} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
3158  fhECharged->SetXTitle("#it{E} (GeV)");
3159  outputContainer->Add(fhECharged);
3160 
3161  fhPtCharged = new TH1F ("hPtCharged","#it{p}_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
3162  fhPtCharged->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3163  outputContainer->Add(fhPtCharged);
3164 
3165  fhPhiCharged = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
3166  fhPhiCharged->SetXTitle("#phi (rad)");
3167  outputContainer->Add(fhPhiCharged);
3168 
3169  fhEtaCharged = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
3170  fhEtaCharged->SetXTitle("#eta ");
3171  outputContainer->Add(fhEtaCharged);
3172 
3173  if(fFillAllTH3)
3174  {
3175  fhEtaPhiECharged = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
3176  netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
3177  fhEtaPhiECharged->SetXTitle("#eta ");
3178  fhEtaPhiECharged->SetYTitle("#phi ");
3179  fhEtaPhiECharged->SetZTitle("#it{E} (GeV) ");
3180  outputContainer->Add(fhEtaPhiECharged);
3181  }
3182  else
3183  {
3184  fhEtaPhiCharged = new TH2F ("hEtaPhiCharged","#eta vs #phi for #it{E} > 0.5 GeV, reconstructed clusters, with matched track",
3185  netabins,etamin,etamax,nphibins,phimin,phimax);
3186  fhEtaPhiCharged->SetXTitle("#eta ");
3187  fhEtaPhiCharged->SetYTitle("#phi (rad)");
3188  outputContainer->Add(fhEtaPhiCharged);
3189  }
3190 
3191  fh1EOverP = new TH2F("h1EOverP","TRACK matches #it{E}/#it{p}",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3192  fh1EOverP->SetYTitle("#it{E}/#it{p}");
3193  fh1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3194  outputContainer->Add(fh1EOverP);
3195 
3196  fh2dR = new TH2F("h2dR","TRACK matches #Delta #it{R}",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
3197  fh2dR->SetYTitle("#Delta #it{R} (rad)");
3198  fh2dR->SetXTitle("#it{E} cluster (GeV)");
3199  outputContainer->Add(fh2dR) ;
3200 
3201  fh2MatchdEdx = new TH2F("h2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
3202  fh2MatchdEdx->SetXTitle("p (GeV/#it{c})");
3203  fh2MatchdEdx->SetYTitle("<#it{dE/dx}>");
3204  outputContainer->Add(fh2MatchdEdx);
3205 
3206  fh2EledEdx = new TH2F("h2EledEdx","#it{dE/dx} vs. #it{p} for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
3207  fh2EledEdx->SetXTitle("p (GeV/#it{c})");
3208  fh2EledEdx->SetYTitle("<#it{dE/dx}>");
3209  outputContainer->Add(fh2EledEdx) ;
3210 
3211  fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches #it{E}/#it{p}, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3212  fh1EOverPR02->SetYTitle("#it{E}/#it{p}");
3213  fh1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3214  outputContainer->Add(fh1EOverPR02);
3215 
3216  fh1EleEOverP = new TH2F("h1EleEOverP","Electron candidates #it{E}/#it{p} (60<#it{dE/dx}<100)",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3217  fh1EleEOverP->SetYTitle("#it{E}/#it{p}");
3218  fh1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3219  outputContainer->Add(fh1EleEOverP);
3220 
3221 
3222  // Projections per SM
3223 
3224  fh1EOverPMod = new TH2F("h1EOverP_PerModule","TRACK matches #it{E}/#it{p}, #it{E}_{cl}&#it{p}_{tr}>0.5 Gev/#it{c}",nPoverEbins,eOverPmin,eOverPmax,fNModules,0,fNModules);
3225  fh1EOverPMod->SetXTitle("#it{E}/#it{p}");
3226  fh1EOverPMod->SetYTitle("Module");
3227  outputContainer->Add(fh1EOverPMod);
3228 
3229  fh2dRMod = new TH2F("h2dR_PerModule","TRACK matches #Delta #it{R}, #it{E}_{cl}&#it{p}_{tr}>0.5 Gev/#it{c}",ndRbins,dRmin,dRmax,fNModules,0,fNModules);
3230  fh2dRMod->SetXTitle("#Delta #it{R} (rad)");
3231  fh2dRMod->SetYTitle("Module");
3232  outputContainer->Add(fh2dRMod) ;
3233 
3234  fh2MatchdEdxMod = new TH2F("h2MatchdEdx_PerModule","#it{dE/dx} vs. #it{p} for all matches, #it{E}_{cl}&#it{p}_{tr}>0.5 Gev/#it{c}",ndedxbins,dedxmin,dedxmax,fNModules,0,fNModules);
3235  fh2MatchdEdxMod->SetYTitle("Module");
3236  fh2MatchdEdxMod->SetXTitle("<#it{dE/dx}>");
3237  outputContainer->Add(fh2MatchdEdxMod);
3238 
3239  fh2EledEdxMod = new TH2F("h2EledEdx_PerModule","#it{dE/dx} vs. #it{p} for electrons, #it{E}_{cl}&#it{p}_{tr}>0.5 Gev/#it{c}",ndedxbins,dedxmin,dedxmax,fNModules,0,fNModules);
3240  fh2EledEdxMod->SetYTitle("Module");
3241  fh2EledEdxMod->SetXTitle("<#it{dE/dx}>");
3242  outputContainer->Add(fh2EledEdxMod) ;
3243 
3244  fh1EOverPR02Mod = new TH2F("h1EOverPR02_PerModule","TRACK matches #it{E}/#it{p}, all, #it{E}_{cl}&#it{p}_{tr}>0.5 Gev/#it{c}",nPoverEbins,eOverPmin,eOverPmax,fNModules,0,fNModules);
3245  fh1EOverPR02Mod->SetXTitle("#it{E}/#it{p}");
3246  fh1EOverPR02Mod->SetYTitle("Module");
3247  outputContainer->Add(fh1EOverPR02Mod);
3248 
3249  fh1EleEOverPMod = new TH2F("h1EleEOverP_PerModule","Electron candidates #it{E}/#it{p} (60<#it{dE/dx}<100), #it{E}_{cl}&#it{p}_{tr}>0.5 Gev/#it{c}",nPoverEbins,eOverPmin,eOverPmax,fNModules,0,fNModules);
3250  fh1EleEOverPMod->SetXTitle("#it{E}/#it{p}");
3251  fh1EleEOverPMod->SetYTitle("Module");
3252  outputContainer->Add(fh1EleEOverPMod);
3253  }
3254 
3255  if(fFillAllPi0Histo)
3256  {
3257  fhIM = new TH2F ("hIM","Cluster pairs (EMCAL or PHOS) Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1",
3258  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3259  fhIM->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3260  fhIM->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3261  outputContainer->Add(fhIM);
3262 
3263  fhIMDiff = new TH2F ("hIMDiff","Cluster pairs (EMCAL or PHOS) Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1, different SM",
3264  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3265  fhIMDiff->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3266  fhIMDiff->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3267  outputContainer->Add(fhIMDiff);
3268 
3269  fhIMSame = new TH2F ("hIMSame","Cluster pairs (EMCAL or PHOS) Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1, same SM",
3270  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3271  fhIMSame->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3272  fhIMSame->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3273  outputContainer->Add(fhIMSame);
3274 
3275  if(fNModules > 12 && (GetCalorimeter() == kEMCAL || GetCalorimeter() == kDCAL))
3276  {
3277  fhIMDCAL = new TH2F ("hIMDCAL","Cluster pairs in DCAL Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1",
3278  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3279  fhIMDCAL->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3280  fhIMDCAL->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3281  outputContainer->Add(fhIMDCAL);
3282 
3283  fhIMDCALDiff = new TH2F ("hIMDCALDiff","Cluster pairs in DCAL Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1, different SM",
3284  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3285  fhIMDCALDiff->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3286  fhIMDCALDiff->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3287  outputContainer->Add(fhIMDCALDiff);
3288 
3289  fhIMDCALSame = new TH2F ("hIMDCALSame","Cluster pairs in DCAL Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1, same SM",
3290  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3291  fhIMDCALSame->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3292  fhIMDCALSame->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3293  outputContainer->Add(fhIMDCALSame);
3294 
3295  fhIMDCALPHOS = new TH2F ("hIMDCALPHOS","Cluster pairs in DCAL-PHOS Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1",
3296  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3297  fhIMDCALPHOS->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3298  fhIMDCALPHOS->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3299  outputContainer->Add(fhIMDCALPHOS);
3300 
3301  fhIMDCALPHOSSame = new TH2F ("hIMDCALPHOSSame","Cluster pairs in DCAL-PHOS Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1, same #phi sector",
3302  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3303  fhIMDCALPHOSSame->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3304  fhIMDCALPHOSSame->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3305  outputContainer->Add(fhIMDCALPHOSSame);
3306  }
3307 
3308  if(fFillInvMassInEMCALWithPHOSDCalAcc)
3309  {
3310  fhIMEMCALPHOS = new TH2F ("hIMEMCALPHOS","Cluster pairs in DCAL-PHOS Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1",
3311  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3312  fhIMEMCALPHOS->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3313  fhIMEMCALPHOS->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3314  outputContainer->Add(fhIMEMCALPHOS);
3315 
3316  fhIMEMCALPHOSSame = new TH2F ("hIMEMCALPHOSSame","Cluster pairs in DCAL-PHOS Invariant mass vs reconstructed pair #it{p}_{T}, ncell > 1, same #phi sector",
3317  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3318  fhIMEMCALPHOSSame->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3319  fhIMEMCALPHOSSame->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3320  outputContainer->Add(fhIMEMCALPHOSSame);
3321  }
3322 
3323  if(fFillPi0PairDiffTime)
3324  {
3325  fhClusterPairDiffTimeEPi0Mass = new TH2F("hClusterPairDiffTimeEPi0Mass","cluster pair time difference vs E, only good clusters",
3326  nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
3327  fhClusterPairDiffTimeEPi0Mass->SetXTitle("#it{E}_{cluster} (GeV)");
3328  fhClusterPairDiffTimeEPi0Mass->SetYTitle("#Delta #it{t} (ns)");
3329  outputContainer->Add(fhClusterPairDiffTimeEPi0Mass);
3330 
3331  fhClusterPairDiffTimeEPi0MassSame = new TH2F("hClusterPairDiffTimeEPi0MassSame","cluster pair time difference vs E, only good clusters",
3332  nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
3333  fhClusterPairDiffTimeEPi0MassSame->SetXTitle("#it{E}_{cluster} (GeV)");
3334  fhClusterPairDiffTimeEPi0MassSame->SetYTitle("#Delta #it{t} (ns)");
3335  outputContainer->Add(fhClusterPairDiffTimeEPi0MassSame);
3336 
3337  if(fNModules > 12 && (GetCalorimeter() == kEMCAL || GetCalorimeter() == kDCAL))
3338  {
3339  fhClusterPairDiffTimeEPi0MassDCal = new TH2F("hClusterPairDiffTimeEPi0MassDCal","cluster pair time difference vs E, only good clusters",
3340  nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
3341  fhClusterPairDiffTimeEPi0MassDCal->SetXTitle("#it{E}_{cluster} (GeV)");
3342  fhClusterPairDiffTimeEPi0MassDCal->SetYTitle("#Delta #it{t} (ns)");
3343  outputContainer->Add(fhClusterPairDiffTimeEPi0MassDCal);
3344 
3345  fhClusterPairDiffTimeEPi0MassDCalSame = new TH2F("hClusterPairDiffTimeEPi0MassDCalSame","cluster pair time difference vs E, only good clusters",
3346  nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
3347  fhClusterPairDiffTimeEPi0MassDCalSame->SetXTitle("#it{E}_{cluster} (GeV)");
3348  fhClusterPairDiffTimeEPi0MassDCalSame->SetYTitle("#Delta #it{t} (ns)");
3349  outputContainer->Add(fhClusterPairDiffTimeEPi0MassDCalSame);
3350  }
3351  }
3352 
3353  fhAsym = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
3354  fhAsym->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3355  fhAsym->SetYTitle("#it{Asymmetry}");
3356  outputContainer->Add(fhAsym);
3357 
3358  if(fFillInvMassOpenAngle)
3359  {
3360  fhOpAngle = new TH2F ("hOpeningAngle","Cluster pairs opening angle vs reconstructed pair #it{p}_{T}, ncell > 1",
3361  nptbins,ptmin,ptmax, 180,0,TMath::Pi());
3362  fhOpAngle->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3363  fhOpAngle->SetYTitle("Opening angle (degrees)");
3364  outputContainer->Add(fhOpAngle);
3365 
3366  Int_t nBinOpAngle = TMath::Nint(fInvMassMaxOpenAngle*TMath::RadToDeg()); // bin of 1 degree size
3367  fhIMvsOpAngle = new TH2F ("hIMvsOpAngle","Cluster pairs Invariant mass vs reconstructed pair opening angle, ncell > 1",
3368  nBinOpAngle,0.,fInvMassMaxOpenAngle,nmassbins,massmin,massmax);
3369  fhIMvsOpAngle->SetXTitle("Opening angle (degrees)");
3370  fhIMvsOpAngle->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
3371  outputContainer->Add(fhIMvsOpAngle);
3372  }
3373  }
3374 
3375  if(fFillAllPosHisto2)
3376  {
3377  fhXYZ = new TH3F ("hXYZ","Cluster: #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
3378  fhXYZ->SetXTitle("#it{x} (cm)");
3379  fhXYZ->SetYTitle("#it{y} (cm)");
3380  fhXYZ->SetZTitle("#it{z} (cm) ");
3381  outputContainer->Add(fhXYZ);
3382 
3383  fhXE = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
3384  fhXE->SetXTitle("#it{x} (cm)");
3385  fhXE->SetYTitle("#it{E} (GeV)");
3386  outputContainer->Add(fhXE);
3387 
3388  fhYE = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
3389  fhYE->SetXTitle("#it{y} (cm)");
3390  fhYE->SetYTitle("#it{E} (GeV)");
3391  outputContainer->Add(fhYE);
3392 
3393  fhZE = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
3394  fhZE->SetXTitle("#it{z} (cm)");
3395  fhZE->SetYTitle("#it{E} (GeV)");
3396  outputContainer->Add(fhZE);
3397 
3398  fhRE = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
3399  fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
3400  fhRE->SetYTitle("#it{E} (GeV)");
3401  outputContainer->Add(fhRE);
3402 
3403  fhXNCells = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
3404  fhXNCells->SetXTitle("#it{x} (cm)");
3405  fhXNCells->SetYTitle("N cells per cluster");
3406  outputContainer->Add(fhXNCells);
3407 
3408  fhYNCells = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
3409  fhYNCells->SetXTitle("#it{y} (cm)");
3410  fhYNCells->SetYTitle("N cells per cluster");
3411  outputContainer->Add(fhYNCells);
3412 
3413  fhZNCells = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
3414  fhZNCells->SetXTitle("#it{z} (cm)");
3415  fhZNCells->SetYTitle("N cells per cluster");
3416  outputContainer->Add(fhZNCells);
3417 
3418  fhRNCells = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
3419  fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
3420  fhRNCells->SetYTitle("N cells per cluster");
3421  outputContainer->Add(fhRNCells);
3422  }
3423 
3424  if(fFillAllPosHisto)
3425  {
3426  fhRCellE = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
3427  fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
3428  fhRCellE->SetYTitle("#it{E} (GeV)");
3429  outputContainer->Add(fhRCellE);
3430 
3431  fhXCellE = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
3432  fhXCellE->SetXTitle("#it{x} (cm)");
3433  fhXCellE->SetYTitle("#it{E} (GeV)");
3434  outputContainer->Add(fhXCellE);
3435 
3436  fhYCellE = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
3437  fhYCellE->SetXTitle("#it{y} (cm)");
3438  fhYCellE->SetYTitle("#it{E} (GeV)");
3439  outputContainer->Add(fhYCellE);
3440 
3441  fhZCellE = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
3442  fhZCellE->SetXTitle("#it{z} (cm)");
3443  fhZCellE->SetYTitle("#it{E} (GeV)");
3444  outputContainer->Add(fhZCellE);
3445 
3446  fhXYZCell = new TH3F ("hXYZCell","Cell : #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
3447  fhXYZCell->SetXTitle("#it{x} (cm)");
3448  fhXYZCell->SetYTitle("#it{y} (cm)");
3449  fhXYZCell->SetZTitle("#it{z} (cm)");
3450  outputContainer->Add(fhXYZCell);
3451 
3452  Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
3453  Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
3454  Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
3455  Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
3456 
3457  fhDeltaCellClusterRNCells = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
3458  fhDeltaCellClusterRNCells->SetXTitle("#it{r} = #sqrt{x^{2}+y^{2}}, #it{r}_{clus}-#it{r}_{cell} (cm)");
3459  fhDeltaCellClusterRNCells->SetYTitle("#it{n}_{cells per cluster}");
3460  outputContainer->Add(fhDeltaCellClusterRNCells);
3461 
3462  fhDeltaCellClusterXNCells = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
3463  fhDeltaCellClusterXNCells->SetXTitle("#it{x}_{clus}-#it{x}_{cell} (cm)");
3464  fhDeltaCellClusterXNCells->SetYTitle("#it{n}_{cells per cluster}");
3465  outputContainer->Add(fhDeltaCellClusterXNCells);
3466 
3467  fhDeltaCellClusterYNCells = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
3468  fhDeltaCellClusterYNCells->SetXTitle("#it{y}_{clus}-#it{y}_{cell} (cm)");
3469  fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
3470  outputContainer->Add(fhDeltaCellClusterYNCells);
3471 
3472  fhDeltaCellClusterZNCells = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
3473  fhDeltaCellClusterZNCells->SetXTitle("#it{z}_{clus}-#it{z}_{cell} (cm)");
3474  fhDeltaCellClusterZNCells->SetYTitle("#it{n}_{cells per cluster}");
3475  outputContainer->Add(fhDeltaCellClusterZNCells);
3476 
3477  fhDeltaCellClusterRE = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
3478  fhDeltaCellClusterRE->SetXTitle("#it{r} = #sqrt{x^{2}+y^{2}}, #it{r}_{clus}-#it{r}_{cell} (cm)");
3479  fhDeltaCellClusterRE->SetYTitle("#it{E} (GeV)");
3480  outputContainer->Add(fhDeltaCellClusterRE);
3481 
3482  fhDeltaCellClusterXE = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
3483  fhDeltaCellClusterXE->SetXTitle("#it{x}_{clus}-#it{x}_{cell} (cm)");
3484  fhDeltaCellClusterXE->SetYTitle("#it{E} (GeV)");
3485  outputContainer->Add(fhDeltaCellClusterXE);
3486 
3487  fhDeltaCellClusterYE = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
3488  fhDeltaCellClusterYE->SetXTitle("#it{y}_{clus}-#it{y}_{cell} (cm)");
3489  fhDeltaCellClusterYE->SetYTitle("#it{E} (GeV)");
3490  outputContainer->Add(fhDeltaCellClusterYE);
3491 
3492  fhDeltaCellClusterZE = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
3493  fhDeltaCellClusterZE->SetXTitle("#it{z}_{clus}-#it{z}_{cell} (cm)");
3494  fhDeltaCellClusterZE->SetYTitle("#it{E} (GeV)");
3495  outputContainer->Add(fhDeltaCellClusterZE);
3496 
3497  if(fFillAllTH3)
3498  {
3499  fhEtaPhiAmpCell = new TH3F ("hEtaPhiAmpCell","Cell #eta vs cell #phi vs cell energy",
3500  netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
3501  fhEtaPhiAmpCell->SetXTitle("#eta ");
3502  fhEtaPhiAmpCell->SetYTitle("#phi (rad)");
3503  fhEtaPhiAmpCell->SetZTitle("#it{E} (GeV) ");
3504  outputContainer->Add(fhEtaPhiAmpCell);
3505  }
3506  else
3507  {
3508  fhEtaPhiCell = new TH2F ("hEtaPhiCell","Cell #eta vs cell #phi vs cell energy",
3509  netabins,etamin,etamax,nphibins,phimin,phimax);
3510  fhEtaPhiCell->SetXTitle("#eta ");
3511  fhEtaPhiCell->SetYTitle("#phi (rad)");
3512  outputContainer->Add(fhEtaPhiCell);
3513  }
3514  }
3515 
3516  // Calorimeter cells
3517 
3518  fhNCells = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
3519  fhNCells->SetXTitle("#it{n}_{cells}");
3520  outputContainer->Add(fhNCells);
3521 
3522  fhNCellsCutAmpMin = new TH1F ("hNCellsCutAmpMin",Form("# cells amp > %1.2f-%1.2f",fEMCALCellAmpMin,fPHOSCellAmpMin), ncebins,ncemin+0.5,ncemax);
3523  fhNCellsCutAmpMin->SetXTitle("#it{n}_{cells}");
3524  outputContainer->Add(fhNCellsCutAmpMin);
3525 
3526  fhAmplitude = new TH1F ("hAmplitude","#it{E}_{cell}", nptbins*2,ptmin,ptmax);
3527  fhAmplitude->SetXTitle("#it{E}_{cell} (GeV)");
3528  outputContainer->Add(fhAmplitude);
3529 
3530  fhAmpId = new TH2F ("hAmpId","#it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
3531  fhAmpId->SetXTitle("#it{E}_{cell} (GeV)");
3532  outputContainer->Add(fhAmpId);
3533 
3534  fhAmpIdLowGain = new TH2F ("hAmpIdLG","Low gain: #it{E}_{cell}", nfineptbins,15,ptfinemax+15,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
3535  fhAmpIdLowGain->SetXTitle("#it{E}_{cell} (GeV)");
3536  outputContainer->Add(fhAmpIdLowGain);
3537 
3538  if(fFillAllCellTimeHisto)
3539  {
3540  fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
3541  fhCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
3542  fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max-i} (ns)");
3543  outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
3544 
3545 // fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
3546 // fhClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
3547 // fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
3548 // outputContainer->Add(fhClusterMaxCellDiffAverageTime);
3549 //
3550 // fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
3551 // fhClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
3552 // fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
3553 // outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
3554 
3555  fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
3556  fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
3557  fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
3558  outputContainer->Add(fhCellIdCellLargeTimeSpread);
3559 
3560  fhTime = new TH1F ("hTime","#it{t}_{cell}",ntimebins,timemin,timemax);
3561  fhTime->SetXTitle("#it{t}_{cell} (ns)");
3562  outputContainer->Add(fhTime);
3563 
3564 // fhTimeVz = new TH2F ("hTimeVz","#it{t}_{cell} vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
3565 // fhTimeVz->SetXTitle("|v_{z}| (cm)");
3566 // fhTimeVz->SetYTitle("#it{t}_{cell} (ns)");
3567 // outputContainer->Add(fhTimeVz);
3568 
3569  fhTimeId = new TH2F ("hTimeId","#it{t}_{cell} vs Absolute Id",
3570  ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
3571  fhTimeId->SetXTitle("#it{t}_{cell} (ns)");
3572  fhTimeId->SetYTitle("Cell Absolute Id");
3573  outputContainer->Add(fhTimeId);
3574 
3575  fhTimeAmp = new TH2F ("hTimeAmp","#it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
3576  fhTimeAmp->SetYTitle("#it{t}_{cell} (ns)");
3577  fhTimeAmp->SetXTitle("#it{E}_{cell} (GeV)");
3578  outputContainer->Add(fhTimeAmp);
3579 
3580  for(Int_t bc = 0; bc < 4; bc++)
3581  {
3582  fhTimePerSMPerBC[bc] = new TH2F (Form("hTimePerSM_BC%d",bc),
3583  Form("#it{t}_{cell} vs super-module, for BC/4=%d",bc),
3584  ntimebins,timemin,timemax,fNModules,0,fNModules);
3585  fhTimePerSMPerBC[bc]->SetXTitle("#it{t}_{cell} (ns)");
3586  fhTimePerSMPerBC[bc]->SetYTitle("Module");
3587  outputContainer->Add(fhTimePerSMPerBC[bc]);
3588  }
3589 
3590  fhTimeIdLowGain = new TH2F ("hTimeIdLG","Low gain: #it{t}_{cell} vs Absolute Id",
3591  ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
3592  fhTimeIdLowGain->SetXTitle("#it{t}_{cell} (ns)");
3593  fhTimeIdLowGain->SetYTitle("Cell Absolute Id");
3594  outputContainer->Add(fhTimeIdLowGain);
3595 
3596  fhTimeAmpLowGain = new TH2F ("hTimeAmpLG","Low gain: #it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
3597  fhTimeAmpLowGain->SetYTitle("#it{t}_{cell} (ns)");
3598  fhTimeAmpLowGain->SetXTitle("#it{E}_{cell} (GeV)");
3599  outputContainer->Add(fhTimeAmpLowGain);
3600  }
3601 
3602  fhCellECross = new TH2F ("hCellECross","1 - Energy in cross around cell / cell energy",
3603  nptbins,ptmin,ptmax, 400,-1,1.);
3604  fhCellECross->SetXTitle("#it{E}_{cell} (GeV) ");
3605  fhCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell}");
3606  outputContainer->Add(fhCellECross);
3607 
3608 
3609  if(fCorrelate)
3610  {
3611  // PHOS vs EMCAL
3612  fhEMCALPHOSCorrNClusters = new TH2F ("hEMCALPHOSCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
3613  fhEMCALPHOSCorrNClusters->SetXTitle("number of clusters in EMCAL");
3614  fhEMCALPHOSCorrNClusters->SetYTitle("number of clusters in PHOS");
3615  outputContainer->Add(fhEMCALPHOSCorrNClusters);
3616 
3617  fhEMCALPHOSCorrEClusters = new TH2F ("hEMCALPHOSCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
3618  fhEMCALPHOSCorrEClusters->SetXTitle("#Sigma #it{E} of clusters in EMCAL (GeV)");
3619  fhEMCALPHOSCorrEClusters->SetYTitle("#Sigma #it{E} of clusters in PHOS (GeV)");
3620  outputContainer->Add(fhEMCALPHOSCorrEClusters);
3621 
3622  fhEMCALPHOSCorrNCells = new TH2F ("hEMCALPHOSCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
3623  fhEMCALPHOSCorrNCells->SetXTitle("number of Cells in EMCAL");
3624  fhEMCALPHOSCorrNCells->SetYTitle("number of Cells in PHOS");
3625  outputContainer->Add(fhEMCALPHOSCorrNCells);
3626 
3627  fhEMCALPHOSCorrECells = new TH2F ("hEMCALPHOSCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
3628  fhEMCALPHOSCorrECells->SetXTitle("#Sigma #it{E} of Cells in EMCAL (GeV)");
3629  fhEMCALPHOSCorrECells->SetYTitle("#Sigma #it{E} of Cells in PHOS (GeV)");
3630  outputContainer->Add(fhEMCALPHOSCorrECells);
3631 
3632  // DCal vs EMCAL
3633  fhEMCALDCALCorrNClusters = new TH2F ("hEMCALDCALCorrNClusters","# clusters in EMCAL vs DCAL", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
3634  fhEMCALDCALCorrNClusters->SetXTitle("number of clusters in EMCAL");
3635  fhEMCALDCALCorrNClusters->SetYTitle("number of clusters in DCAL");
3636  outputContainer->Add(fhEMCALDCALCorrNClusters);
3637 
3638  fhEMCALDCALCorrEClusters = new TH2F ("hEMCALDCALCorrEClusters","summed energy of clusters in EMCAL vs DCAL", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
3639  fhEMCALDCALCorrEClusters->SetXTitle("#Sigma #it{E} of clusters in EMCAL (GeV)");
3640  fhEMCALDCALCorrEClusters->SetYTitle("#Sigma #it{E} of clusters in DCAL (GeV)");
3641  outputContainer->Add(fhEMCALDCALCorrEClusters);
3642 
3643  fhEMCALDCALCorrNCells = new TH2F ("hEMCALDCALCorrNCells","# Cells in EMCAL vs DCAL", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
3644  fhEMCALDCALCorrNCells->SetXTitle("number of Cells in EMCAL");
3645  fhEMCALDCALCorrNCells->SetYTitle("number of Cells in DCAL");
3646  outputContainer->Add(fhEMCALDCALCorrNCells);
3647 
3648  fhEMCALDCALCorrECells = new TH2F ("hEMCALDCALCorrECells","summed energy of Cells in EMCAL vs DCAL", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
3649  fhEMCALDCALCorrECells->SetXTitle("#Sigma #it{E} of Cells in EMCAL (GeV)");
3650  fhEMCALDCALCorrECells->SetYTitle("#Sigma #it{E} of Cells in DCAL (GeV)");
3651  outputContainer->Add(fhEMCALDCALCorrECells);
3652 
3653 
3654  // DCAL vs PHOS
3655  fhDCALPHOSCorrNClusters = new TH2F ("hDCALPHOSCorrNClusters","# clusters in DCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
3656  fhDCALPHOSCorrNClusters->SetXTitle("number of clusters in DCAL");
3657  fhDCALPHOSCorrNClusters->SetYTitle("number of clusters in PHOS");
3658  outputContainer->Add(fhDCALPHOSCorrNClusters);
3659 
3660  fhDCALPHOSCorrEClusters = new TH2F ("hDCALPHOSCorrEClusters","summed energy of clusters in DCAL vs PHOS", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
3661  fhDCALPHOSCorrEClusters->SetXTitle("#Sigma #it{E} of clusters in DCAL (GeV)");
3662  fhDCALPHOSCorrEClusters->SetYTitle("#Sigma #it{E} of clusters in PHOS (GeV)");
3663  outputContainer->Add(fhDCALPHOSCorrEClusters);
3664 
3665  fhDCALPHOSCorrNCells = new TH2F ("hDCALPHOSCorrNCells","# Cells in DCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
3666  fhDCALPHOSCorrNCells->SetXTitle("number of Cells in DCAL");
3667  fhDCALPHOSCorrNCells->SetYTitle("number of Cells in PHOS");
3668  outputContainer->Add(fhDCALPHOSCorrNCells);
3669 
3670  fhDCALPHOSCorrECells = new TH2F ("hDCALPHOSCorrECells","summed energy of Cells in DCAL vs PHOS", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
3671  fhDCALPHOSCorrECells->SetXTitle("#Sigma #it{E} of Cells in DCAL (GeV)");
3672  fhDCALPHOSCorrECells->SetYTitle("#Sigma #it{E} of Cells in PHOS (GeV)");
3673  outputContainer->Add(fhDCALPHOSCorrECells);
3674 
3675  // Calorimeter vs. V0 signal
3676 
3677  fhCaloV0SCorrNClusters = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
3678  fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
3679  fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
3680  outputContainer->Add(fhCaloV0SCorrNClusters);
3681 
3682  fhCaloV0SCorrEClusters = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
3683  fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
3684  fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
3685  outputContainer->Add(fhCaloV0SCorrEClusters);
3686 
3687  fhCaloV0SCorrNCells = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
3688  fhCaloV0SCorrNCells->SetXTitle("V0 signal");
3689  fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
3690  outputContainer->Add(fhCaloV0SCorrNCells);
3691 
3692  fhCaloV0SCorrECells = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
3693  fhCaloV0SCorrECells->SetXTitle("V0 signal");
3694  fhCaloV0SCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
3695  outputContainer->Add(fhCaloV0SCorrECells);
3696 
3697  // Calorimeter vs V0 multiplicity
3698 
3699  fhCaloV0MCorrNClusters = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
3700  fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
3701  fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
3702  outputContainer->Add(fhCaloV0MCorrNClusters);
3703 
3704  fhCaloV0MCorrEClusters = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
3705  fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
3706  fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
3707  outputContainer->Add(fhCaloV0MCorrEClusters);
3708 
3709  fhCaloV0MCorrNCells = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
3710  fhCaloV0MCorrNCells->SetXTitle("V0 signal");
3711  fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
3712  outputContainer->Add(fhCaloV0MCorrNCells);
3713 
3714  fhCaloV0MCorrECells = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
3715  fhCaloV0MCorrECells->SetXTitle("V0 signal");
3716  fhCaloV0MCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
3717  outputContainer->Add(fhCaloV0MCorrECells);
3718 
3719  //Calorimeter VS Track multiplicity
3720  fhCaloTrackMCorrNClusters = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
3721  fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
3722  fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
3723  outputContainer->Add(fhCaloTrackMCorrNClusters);
3724 
3725  fhCaloTrackMCorrEClusters = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
3726  fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
3727  fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
3728  outputContainer->Add(fhCaloTrackMCorrEClusters);
3729 
3730  fhCaloTrackMCorrNCells = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
3731  fhCaloTrackMCorrNCells->SetXTitle("# tracks");
3732  fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
3733  outputContainer->Add(fhCaloTrackMCorrNCells);
3734 
3735  fhCaloTrackMCorrECells = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
3736  fhCaloTrackMCorrECells->SetXTitle("# tracks");
3737  fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
3738  outputContainer->Add(fhCaloTrackMCorrECells);
3739 
3740  fhCaloCenNClusters = new TH2F ("hCaloCenNClusters","# clusters in calorimeter vs centrality",100,0,100,nclbins,nclmin,nclmax);
3741  fhCaloCenNClusters->SetYTitle("number of clusters in calorimeter");
3742  fhCaloCenNClusters->SetXTitle("Centrality");
3743  outputContainer->Add(fhCaloCenNClusters);
3744 
3745  fhCaloCenEClusters = new TH2F ("hCaloCenEClusters","summed energy of clusters in calorimeter vs centrality",100,0,100,nptbins,ptmin,ptmax*2);
3746  fhCaloCenEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
3747  fhCaloCenEClusters->SetXTitle("Centrality");
3748  outputContainer->Add(fhCaloCenEClusters);
3749 
3750  fhCaloCenNCells = new TH2F ("hCaloCenNCells","# Cells in calorimeter vs centrality",100,0,100,ncebins,ncemin,ncemax);
3751  fhCaloCenNCells->SetYTitle("number of Cells in calorimeter");
3752  fhCaloCenNCells->SetXTitle("Centrality");
3753  outputContainer->Add(fhCaloCenNCells);
3754 
3755  fhCaloCenECells = new TH2F ("hCaloCenECells","summed energy of Cells in calorimeter vs centrality",100,0,100,nptbins*2,ptmin,ptmax*4);
3756  fhCaloCenECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
3757  fhCaloCenECells->SetXTitle("Centrality");
3758  outputContainer->Add(fhCaloCenECells);
3759 
3760  fhCaloEvPNClusters = new TH2F ("hCaloEvPNClusters","# clusters in calorimeter vs event plane angle",100,0,TMath::Pi(),nclbins,nclmin,nclmax);
3761  fhCaloEvPNClusters->SetYTitle("number of clusters in calorimeter");
3762  fhCaloEvPNClusters->SetXTitle("Event plane angle (rad)");
3763  outputContainer->Add(fhCaloEvPNClusters);
3764 
3765  fhCaloEvPEClusters = new TH2F ("hCaloEvPEClusters","summed energy of clusters in calorimeter vs event plane angle",100,0,TMath::Pi(),nptbins,ptmin,ptmax*2);
3766  fhCaloEvPEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
3767  fhCaloEvPEClusters->SetXTitle("Event plane angle (rad)");
3768  outputContainer->Add(fhCaloEvPEClusters);
3769 
3770  fhCaloEvPNCells = new TH2F ("hCaloEvPNCells","# Cells in calorimeter vs event plane angle",100,0,TMath::Pi(),ncebins,ncemin,ncemax);
3771  fhCaloEvPNCells->SetYTitle("number of Cells in calorimeter");
3772  fhCaloEvPNCells->SetXTitle("Event plane angle (rad)");
3773  outputContainer->Add(fhCaloEvPNCells);
3774 
3775  fhCaloEvPECells = new TH2F ("hCaloEvPECells","summed energy of Cells in calorimeter vs event plane angle",100,0,TMath::Pi(),nptbins*2,ptmin,ptmax*4);
3776  fhCaloEvPECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
3777  fhCaloEvPECells->SetXTitle("Event plane angle (rad)");
3778  outputContainer->Add(fhCaloEvPECells);
3779  } // correlate calorimeters
3780 
3781  // Module histograms
3782 
3783  fhEMod = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
3784  fhEMod->SetXTitle("#it{E} (GeV)");
3785  fhEMod->SetYTitle("Module");
3786  outputContainer->Add(fhEMod);
3787 
3788  fhAmpMod = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
3789  fhAmpMod->SetXTitle("#it{E} (GeV)");
3790  fhAmpMod->SetYTitle("Module");
3791  outputContainer->Add(fhAmpMod);
3792 
3793  fhEWeirdMod = new TH2F ("hEWeird_Mod","Cluster reconstructed Energy in each present Module, ridiculously large E",200,0,10000,fNModules,0,fNModules);
3794  fhEWeirdMod->SetXTitle("#it{E} (GeV)");
3795  fhEWeirdMod->SetYTitle("Module");
3796  outputContainer->Add(fhEWeirdMod);
3797 
3798  fhAmpWeirdMod = new TH2F ("hAmpWeird_Mod","Cell energy in each present Module, ridiculously large E",200,0,10000,fNModules,0,fNModules);
3799  fhAmpWeirdMod->SetXTitle("#it{E} (GeV)");
3800  fhAmpWeirdMod->SetYTitle("Module");
3801  outputContainer->Add(fhAmpWeirdMod);
3802 
3803  if(fFillAllCellTimeHisto)
3804  {
3805  fhTimeMod = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
3806  fhTimeMod->SetXTitle("t (ns)");
3807  fhTimeMod->SetYTitle("Module");
3808  outputContainer->Add(fhTimeMod);
3809  }
3810 
3811  fhNClustersMod = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
3812  fhNClustersMod->SetXTitle("number of clusters");
3813  fhNClustersMod->SetYTitle("Module");
3814  outputContainer->Add(fhNClustersMod);
3815 
3816  fhNCellsMod = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
3817  fhNCellsMod->SetXTitle("#it{n}_{cells}");
3818  fhNCellsMod->SetYTitle("Module");
3819  outputContainer->Add(fhNCellsMod);
3820 
3821  fhSumClustersEnergyMod = new TH2F ("hSumClustersEnergy_Mod","# clusters vs Module", 1000, 0, 2000,fNModules,0,fNModules);
3822  fhSumClustersEnergyMod->SetXTitle("#Sigma_{clusters} #it{E} (GeV)");
3823  fhSumClustersEnergyMod->SetYTitle("Module");
3824  outputContainer->Add(fhSumClustersEnergyMod);
3825 
3826  fhSumCellsAmpMod = new TH2F ("hSumCellsAmp_Mod","# cells vs Module", 1000, 0, 2000,fNModules,0,fNModules);
3827  fhSumCellsAmpMod->SetXTitle("#Sigma_{cells} #it{Amp} (GeV)");
3828  fhSumCellsAmpMod->SetYTitle("Module");
3829  outputContainer->Add(fhSumCellsAmpMod);
3830 
3831  Int_t colmaxs = fNMaxCols;
3832  Int_t rowmaxs = fNMaxRows;
3833  if(GetCalorimeter()==kEMCAL)
3834  {
3835  colmaxs=2*fNMaxCols;
3836  rowmaxs=Int_t(fNModules/2)*fNMaxRows;
3837  }
3838  else
3839  {
3840  rowmaxs=fNModules*fNMaxRows;
3841  }
3842 
3843  fhGridCells = new TH2F ("hGridCells",Form("Entries in grid of cells"),
3844  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
3845  fhGridCells->SetYTitle("row (phi direction)");
3846  fhGridCells->SetXTitle("column (eta direction)");
3847  outputContainer->Add(fhGridCells);
3848 
3849  fhGridCellsE = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
3850  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
3851  fhGridCellsE->SetYTitle("row (phi direction)");
3852  fhGridCellsE->SetXTitle("column (eta direction)");
3853  outputContainer->Add(fhGridCellsE);
3854 
3855  fhGridCellsLowGain = new TH2F ("hGridCellsLG",Form("Low gain: Entries in grid of cells"),
3856  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
3857  fhGridCellsLowGain->SetYTitle("row (phi direction)");
3858  fhGridCellsLowGain->SetXTitle("column (eta direction)");
3859  outputContainer->Add(fhGridCellsLowGain);
3860 
3861  fhGridCellsELowGain = new TH2F ("hGridCellsELG","Low gain: Accumulated energy in grid of cells",
3862  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
3863  fhGridCellsELowGain->SetYTitle("row (phi direction)");
3864  fhGridCellsELowGain->SetXTitle("column (eta direction)");
3865  outputContainer->Add(fhGridCellsELowGain);
3866 
3867  if(fFillAllCellTimeHisto)
3868  {
3869  fhGridCellsTime = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
3870  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
3871  fhGridCellsTime->SetYTitle("row (phi direction)");
3872  fhGridCellsTime->SetXTitle("column (eta direction)");
3873  outputContainer->Add(fhGridCellsTime);
3874 
3875  fhGridCellsTimeLowGain = new TH2F ("hGridCellsTimeLG","Low gain: Accumulated time in grid of cells",
3876  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
3877  fhGridCellsTimeLowGain->SetYTitle("row (phi direction)");
3878  fhGridCellsTimeLowGain->SetXTitle("column (eta direction)");
3879  outputContainer->Add(fhGridCellsTimeLowGain);
3880  }
3881 
3882  fhNCellsSumAmpPerMod = new TH2F*[fNModules];
3883 
3884  fhNClustersSumEnergyPerMod = new TH2F*[fNModules];
3885 
3886  fhIMMod = new TH2F*[fNModules];
3887 
3888  fhNCellsPerClusterMod = new TH2F*[fNModules];
3889 
3890  if(fStudyBadClusters)
3891  fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
3892 
3893  if(fFillAllCellTimeHisto)
3894  fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
3895 
3896  for(Int_t imod = 0; imod < fNModules; imod++)
3897  {
3898  fhNCellsSumAmpPerMod[imod] = new TH2F (Form("hNCellsSumAmp_Mod%d",imod),
3899  Form("# cells in SM vs sum of cells energy in Module %d",imod),
3900  nptbins,ptmin,ptmax*4, ncebins,ncemin,ncemax);
3901  fhNCellsSumAmpPerMod[imod]->SetXTitle("#Sigma #it{Amplitude} (GeV)");
3902  fhNCellsSumAmpPerMod[imod]->SetYTitle("#Sigma #it{n}_{cells}");
3903  outputContainer->Add(fhNCellsSumAmpPerMod[imod]);
3904 
3905  fhNClustersSumEnergyPerMod[imod] = new TH2F (Form("hNClustersSumEnergy_Mod%d",imod),
3906  Form("# clusters in SM vs sum of clusters energy in Module %d",imod),
3907  nptbins,ptmin,ptmax*4, nclbins,nclmin,nclmax);
3908  fhNClustersSumEnergyPerMod[imod]->SetXTitle("#Sigma #it{E} (GeV)");
3909  fhNClustersSumEnergyPerMod[imod]->SetYTitle("#Sigma #it{n}_{clusters}");
3910  outputContainer->Add(fhNClustersSumEnergyPerMod[imod]);
3911 
3912  fhNCellsPerClusterMod[imod] = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
3913  Form("# cells per cluster vs cluster energy in Module %d",imod),
3914  nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
3915  fhNCellsPerClusterMod[imod]->SetXTitle("#it{E} (GeV)");
3916  fhNCellsPerClusterMod[imod]->SetYTitle("#it{n}_{cells per cluster}");
3917  outputContainer->Add(fhNCellsPerClusterMod[imod]);
3918 
3919  if(fStudyBadClusters)
3920  {
3921  fhNCellsPerClusterModNoCut[imod] = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
3922  Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
3923  nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
3924  fhNCellsPerClusterModNoCut[imod]->SetXTitle("#it{E} (GeV)");
3925  fhNCellsPerClusterModNoCut[imod]->SetYTitle("#it{n}_{cells per cluster}");
3926  outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
3927  }
3928 
3929  if(fFillAllCellTimeHisto)
3930  {
3931  for(Int_t ircu = 0; ircu < fNRCU; ircu++)
3932  {
3933  fhTimeAmpPerRCU[imod*fNRCU+ircu] = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
3934  Form("#it{E}_{cell} vs #it{t}_{cell} in Module %d, RCU %d ",imod,ircu),
3935  nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
3936  fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("#it{E} (GeV)");
3937  fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("#it{t} (ns)");
3938  outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
3939  }
3940  }
3941 
3942  if(fFillAllPi0Histo)
3943  {
3944  fhIMMod[imod] = new TH2F (Form("hIM_Mod%d",imod),
3945  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
3946  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
3947  fhIMMod[imod]->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
3948  fhIMMod[imod]->SetYTitle("#it{M}_{cluster pairs} (GeV/#it{c}^{2})");
3949  outputContainer->Add(fhIMMod[imod]);
3950  }
3951  }
3952 
3953  // Monte Carlo Histograms
3954 
3955  TString particleName[] = {
3956  "Photon", "Pi0", "Eta",
3957  "Electron", "PhotonConv",
3958  "NeutralHadron", "ChargedHadron" };
3959 
3960  if(IsDataMC())
3961  {
3962  for(Int_t iPart = 0; iPart < 7; iPart++)
3963  {
3964  for(Int_t iCh = 0; iCh < 2; iCh++)
3965  {
3966  fhRecoMCRatioE[iPart][iCh] = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
3967  Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh),
3968  nptbins, ptmin, ptmax, 200,0,2);
3969  fhRecoMCRatioE[iPart][iCh]->SetYTitle("#it{E}_{reconstructed}/#it{E}_{generated}");
3970  fhRecoMCRatioE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
3971  outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
3972 
3973 
3974  fhRecoMCDeltaE[iPart][iCh] = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
3975  Form("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh),
3976  nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
3977  fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta #it{E} (GeV)");
3978  fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
3979  outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
3980 
3981  fhRecoMCDeltaPhi[iPart][iCh] = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
3982  Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
3983  nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
3984  fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
3985  fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
3986  outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
3987 
3988  fhRecoMCDeltaEta[iPart][iCh] = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
3989  Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
3990  nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
3991  fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
3992  fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
3993  outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
3994 
3995  fhRecoMCE[iPart][iCh] = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
3996  Form("#it{E} distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
3997  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3998  fhRecoMCE[iPart][iCh]->SetXTitle("#it{E}_{rec} (GeV)");
3999  fhRecoMCE[iPart][iCh]->SetYTitle("#it{E}_{gen} (GeV)");
4000  outputContainer->Add(fhRecoMCE[iPart][iCh]);
4001 
4002  fhRecoMCPhi[iPart][iCh] = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
4003  Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
4004  nphibins,phimin,phimax, nphibins,phimin,phimax);
4005  fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
4006  fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
4007  outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
4008 
4009  fhRecoMCEta[iPart][iCh] = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
4010  Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
4011  netabins,etamin,etamax,netabins,etamin,etamax);
4012  fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
4013  fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
4014  outputContainer->Add(fhRecoMCEta[iPart][iCh]);
4015  }
4016  }
4017 
4018  // Pure MC
4019 
4020  for(Int_t iPart = 0; iPart < 4; iPart++)
4021  {
4022  fhGenMCE [iPart] = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
4023  Form("#it{E} of generated %s",particleName[iPart].Data()),
4024  nptbins,ptmin,ptmax);
4025 
4026  fhGenMCPt[iPart] = new TH1F(Form("hGenMCPt_%s",particleName[iPart].Data()) ,
4027  Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
4028  nptbins,ptmin,ptmax);
4029 
4030  fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
4031  Form("Y vs #phi of generated %s",particleName[iPart].Data()),
4032  200,-1,1,360,0,TMath::TwoPi());
4033 
4034  fhGenMCE [iPart] ->SetXTitle("#it{E} (GeV)");
4035  fhGenMCPt[iPart] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4036  fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
4037  fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
4038 
4039  outputContainer->Add(fhGenMCE [iPart]);
4040  outputContainer->Add(fhGenMCPt [iPart]);
4041  outputContainer->Add(fhGenMCEtaPhi[iPart]);
4042 
4043 
4044  fhGenMCAccE [iPart] = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
4045  Form("#it{E} of generated %s",particleName[iPart].Data()),
4046  nptbins,ptmin,ptmax);
4047  fhGenMCAccPt[iPart] = new TH1F(Form("hGenMCAccPt_%s",particleName[iPart].Data()) ,
4048  Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
4049  nptbins,ptmin,ptmax);
4050  fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
4051  Form("Y vs #phi of generated %s",particleName[iPart].Data()),
4052  netabins,etamin,etamax,nphibins,phimin,phimax);
4053 
4054  fhGenMCAccE [iPart] ->SetXTitle("#it{E} (GeV)");
4055  fhGenMCAccPt[iPart] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4056  fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
4057  fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
4058 
4059  outputContainer->Add(fhGenMCAccE [iPart]);
4060  outputContainer->Add(fhGenMCAccPt [iPart]);
4061  outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
4062 
4063  }
4064 
4065  // Vertex of generated particles
4066 
4067  fhEMVxyz = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
4068  fhEMVxyz->SetXTitle("#it{v}_{x}");
4069  fhEMVxyz->SetYTitle("#it{v}_{y}");
4070  //fhEMVxyz->SetZTitle("v_{z}");
4071  outputContainer->Add(fhEMVxyz);
4072 
4073  fhHaVxyz = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
4074  fhHaVxyz->SetXTitle("#it{v}_{x}");
4075  fhHaVxyz->SetYTitle("#it{v}_{y}");
4076  //fhHaVxyz->SetZTitle("v_{z}");
4077  outputContainer->Add(fhHaVxyz);
4078 
4079  fhEMR = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
4080  fhEMR->SetXTitle("#it{E} (GeV)");
4081  fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
4082  outputContainer->Add(fhEMR);
4083 
4084  fhHaR = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
4085  fhHaR->SetXTitle("#it{E} (GeV)");
4086  fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
4087  outputContainer->Add(fhHaR);
4088 
4089  // Track Matching
4090 
4091  fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4092  fhMCEle1EOverP->SetYTitle("#it{E}/#it{p}");
4093  fhMCEle1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4094  outputContainer->Add(fhMCEle1EOverP);
4095 
4096  fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
4097  fhMCEle1dR->SetXTitle("#Delta #it{R} (rad)");
4098  outputContainer->Add(fhMCEle1dR) ;
4099 
4100  fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
4101  fhMCEle2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
4102  fhMCEle2MatchdEdx->SetYTitle("<#it{dE/dx}>");
4103  outputContainer->Add(fhMCEle2MatchdEdx);
4104 
4105  fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches #it{E}/#it{p}, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4106  fhMCChHad1EOverP->SetYTitle("#it{E}/#it{p}");
4107  fhMCChHad1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4108  outputContainer->Add(fhMCChHad1EOverP);
4109 
4110  fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
4111  fhMCChHad1dR->SetXTitle("#Delta R (rad)");
4112  outputContainer->Add(fhMCChHad1dR) ;
4113 
4114  fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
4115  fhMCChHad2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
4116  fhMCChHad2MatchdEdx->SetYTitle("#it{dE/dx}>");
4117  outputContainer->Add(fhMCChHad2MatchdEdx);
4118 
4119  fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches #it{E}/#it{p}, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4120  fhMCNeutral1EOverP->SetYTitle("#it{E}/#it{p}");
4121  fhMCNeutral1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4122  outputContainer->Add(fhMCNeutral1EOverP);
4123 
4124  fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
4125  fhMCNeutral1dR->SetXTitle("#Delta #it{R} (rad)");
4126  outputContainer->Add(fhMCNeutral1dR) ;
4127 
4128  fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
4129  fhMCNeutral2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
4130  fhMCNeutral2MatchdEdx->SetYTitle("#it{dE/dx}>");
4131  outputContainer->Add(fhMCNeutral2MatchdEdx);
4132 
4133  fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4134  fhMCEle1EOverPR02->SetYTitle("#it{E}/#it{p}");
4135  fhMCEle1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4136  outputContainer->Add(fhMCEle1EOverPR02);
4137 
4138  fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches #it{E}/#it{p}, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4139  fhMCChHad1EOverPR02->SetYTitle("#it{E}/#it{p}");
4140  fhMCChHad1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4141  outputContainer->Add(fhMCChHad1EOverPR02);
4142 
4143  fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches #it{E}/#it{p}, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4144  fhMCNeutral1EOverPR02->SetYTitle("#it{E}/#it{p}");
4145  fhMCNeutral1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4146  outputContainer->Add(fhMCNeutral1EOverPR02);
4147 
4148  fhMCEle1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4149  fhMCEle1EleEOverP->SetYTitle("#it{E}/#it{p}");
4150  fhMCEle1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4151  outputContainer->Add(fhMCEle1EleEOverP);
4152 
4153  fhMCChHad1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4154  fhMCChHad1EleEOverP->SetYTitle("#it{E}/#it{p}");
4155  fhMCChHad1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4156  outputContainer->Add(fhMCChHad1EleEOverP);
4157 
4158  fhMCNeutral1EleEOverP = new TH2F("hMCNeutral1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
4159  fhMCNeutral1EleEOverP->SetYTitle("#it{E}/#it{p}");
4160  fhMCNeutral1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4161  outputContainer->Add(fhMCNeutral1EleEOverP);
4162  }
4163 
4164 
4165  if(fFillEBinAcceptanceHisto)
4166  {
4167  for(Int_t ie=0; ie<fNEBinCuts; ie++)
4168  {
4169  fhEBinClusterEtaPhi[ie] = new TH2F
4170  (Form("hEBin%d_Cluster_EtaPhi",ie),
4171  Form("#eta vs #phi, cluster, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
4172  netabins,etamin,etamax,nphibins,phimin,phimax);
4173  fhEBinClusterEtaPhi[ie]->SetYTitle("#phi (rad)");
4174  fhEBinClusterEtaPhi[ie]->SetXTitle("#eta");
4175  outputContainer->Add(fhEBinClusterEtaPhi[ie]) ;
4176 
4177  fhEBinClusterColRow[ie] = new TH2F
4178  (Form("hEBin%d_Cluster_ColRow",ie),
4179  Form("column vs row, cluster max E cell, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
4180  96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
4181  fhEBinClusterColRow[ie]->SetYTitle("row");
4182  fhEBinClusterColRow[ie]->SetXTitle("column");
4183  outputContainer->Add(fhEBinClusterColRow[ie]) ;
4184 
4185  fhEBinCellColRow[ie] = new TH2F
4186  (Form("hEBin%d_Cell_ColRow",ie),
4187  Form("column vs row, cell, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fEBinCuts[ie],fEBinCuts[ie+1]),
4188  96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
4189  fhEBinCellColRow[ie]->SetYTitle("row");
4190  fhEBinCellColRow[ie]->SetXTitle("column");
4191  outputContainer->Add(fhEBinCellColRow[ie]) ;
4192  }
4193  }
4194 
4195  if(fStudyM02Dependence)
4196  {
4197  fhClusterTimeEnergyM02 = new TH3F ("hClusterTimeEnergyM02","energy vs TOF vs M02, reconstructed clusters",
4198  11,5.5,16.5,45,-25,20,100,0.,1);
4199  fhClusterTimeEnergyM02->SetXTitle("#it{E} (GeV) ");
4200  fhClusterTimeEnergyM02->SetYTitle("TOF (ns)");
4201  fhClusterTimeEnergyM02->SetZTitle("#lambda_{0}^{2}");
4202  outputContainer->Add(fhClusterTimeEnergyM02);
4203 
4204  if(fFillAllCellTimeHisto)
4205  {
4206  fhCellTimeSpreadRespectToCellMaxM02 = new TH3F ("hCellTimeSpreadRespectToCellMaxM02","t_{cell max}-t_{cell i} vs E vs M02",
4207  11,5.5,16.5,100,-100,100,100,0.,1);
4208  fhCellTimeSpreadRespectToCellMaxM02->SetXTitle("#it{E} (GeV)");
4209  fhCellTimeSpreadRespectToCellMaxM02->SetYTitle("#Delta #it{t}_{cell max-i} (ns)");
4210  fhCellTimeSpreadRespectToCellMaxM02->SetZTitle("#lambda_{0}^{2}");
4211  outputContainer->Add(fhCellTimeSpreadRespectToCellMaxM02);
4212  }
4213 
4214  fhClusterMaxCellCloseCellRatioM02 = new TH3F ("hClusterMaxCellCloseCellRatioM02","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
4215  11,5.5,16.5, 100,0,1.,100,0.,1);
4216  fhClusterMaxCellCloseCellRatioM02->SetXTitle("#it{E}_{cluster} (GeV) ");
4217  fhClusterMaxCellCloseCellRatioM02->SetYTitle("#it{E}_{cell i}/#it{E}_{cell max}");
4218  fhClusterMaxCellCloseCellRatioM02->SetZTitle("#lambda_{0}^{2}");
4219  outputContainer->Add(fhClusterMaxCellCloseCellRatioM02);
4220 
4221  fhClusterMaxCellCloseCellDiffM02 = new TH3F ("hClusterMaxCellCloseCellDiffM02","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
4222  11,5.5,16.5, 100,0,100,100,0.,1);
4223  fhClusterMaxCellCloseCellDiffM02->SetXTitle("#it{E}_{cluster} (GeV) ");
4224  fhClusterMaxCellCloseCellDiffM02->SetYTitle("#it{E}_{cell max}-#it{E}_{cell i} (GeV)");
4225  fhClusterMaxCellCloseCellDiffM02->SetZTitle("#lambda_{0}^{2}");
4226  outputContainer->Add(fhClusterMaxCellCloseCellDiffM02);
4227 
4228  fhClusterMaxCellDiffM02 = new TH3F ("hClusterMaxCellDiffM02","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
4229  11,5.5,16.5, 100,0,1.,100,0.,1);
4230  fhClusterMaxCellDiffM02->SetXTitle("#it{E}_{cluster} (GeV) ");
4231  fhClusterMaxCellDiffM02->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
4232  fhClusterMaxCellDiffM02->SetZTitle("#lambda_{0}^{2}");
4233  outputContainer->Add(fhClusterMaxCellDiffM02);
4234 
4235  fhClusterMaxCellECrossM02 = new TH3F ("hClusterMaxCellECrossM02",
4236  "1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
4237  11,5.5,16.5, 100,0.0,1.,100,0.,1);
4238  fhClusterMaxCellECrossM02->SetXTitle("#it{E}_{cluster} (GeV) ");
4239  fhClusterMaxCellECrossM02->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
4240  fhClusterMaxCellECrossM02->SetZTitle("#lambda_{0}^{2}");
4241  outputContainer->Add(fhClusterMaxCellECrossM02);
4242 
4243  fhNCellsPerClusterM02 = new TH3F ("hNCellsPerClusterM02","# cells per cluster vs energy",
4244  11,5.5,16.5, 35,0.5,35.5,100,0.,1);
4245  fhNCellsPerClusterM02->SetXTitle("#it{E} (GeV)");
4246  fhNCellsPerClusterM02->SetYTitle("#it{n}_{cells}");
4247  fhNCellsPerClusterM02->SetZTitle("#lambda_{0}^{2}");
4248  outputContainer->Add(fhNCellsPerClusterM02);
4249  }
4250 
4251  // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
4252  // printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
4253 
4254  return outputContainer;
4255 }
4256 
4257 //______________________________________________________________________________________
4260 //______________________________________________________________________________________
4261 Float_t AliAnaCalorimeterQA::GetECross(Int_t absID, AliVCaloCells* cells, Float_t dtcut)
4262 {
4263  Int_t icol =-1, irow=-1,iRCU = -1;
4264  Int_t imod = GetModuleNumberCellIndexes(absID, GetCalorimeter(), icol, irow, iRCU);
4265 
4266  if ( GetCalorimeter() == kEMCAL)
4267  {
4268  // Get close cells index, energy and time, not in corners
4269 
4270  Int_t absID1 = -1;
4271  Int_t absID2 = -1;
4272 
4273  if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
4274  if( irow > 0 ) absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
4275 
4276  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
4277  Int_t absID3 = -1;
4278  Int_t absID4 = -1;
4279 
4280  if ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
4281  {
4282  absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
4283  absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol-1);
4284  }
4285  else if( icol == 0 && imod%2 )
4286  {
4287  absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol+1);
4288  absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
4289  }
4290  else
4291  {
4292  if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
4293  absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
4294  if( icol > 0 )
4295  absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
4296  }
4297 
4298  // Recalibrate cell energy if needed
4299  //Float_t ecell = cells->GetCellAmplitude(absID);
4300  //GetCaloUtils()->RecalibrateCellAmplitude(ecell,GetCalorimeter(), absID);
4301  Double_t tcell = cells->GetCellTime(absID);
4302  GetCaloUtils()->RecalibrateCellTime(tcell, GetCalorimeter(), absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
4303 
4304  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
4305  Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
4306 
4307  if(absID1 > 0 )
4308  {
4309  ecell1 = cells->GetCellAmplitude(absID1);
4310  GetCaloUtils()->RecalibrateCellAmplitude(ecell1, GetCalorimeter(), absID1);
4311  tcell1 = cells->GetCellTime(absID1);
4312  GetCaloUtils()->RecalibrateCellTime (tcell1, GetCalorimeter(), absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
4313  }
4314 
4315  if(absID2 > 0 )
4316  {
4317  ecell2 = cells->GetCellAmplitude(absID2);
4318  GetCaloUtils()->RecalibrateCellAmplitude(ecell2, GetCalorimeter(), absID2);
4319  tcell2 = cells->GetCellTime(absID2);
4320  GetCaloUtils()->RecalibrateCellTime (tcell2, GetCalorimeter(), absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
4321  }
4322 
4323  if(absID3 > 0 )
4324  {
4325  ecell3 = cells->GetCellAmplitude(absID3);
4326  GetCaloUtils()->RecalibrateCellAmplitude(ecell3, GetCalorimeter(), absID3);
4327  tcell3 = cells->GetCellTime(absID3);
4328  GetCaloUtils()->RecalibrateCellTime (tcell3, GetCalorimeter(), absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
4329  }
4330 
4331  if(absID4 > 0 )
4332  {
4333  ecell4 = cells->GetCellAmplitude(absID4);
4334  GetCaloUtils()->RecalibrateCellAmplitude(ecell4, GetCalorimeter(), absID4);
4335  tcell4 = cells->GetCellTime(absID4);
4336  GetCaloUtils()->RecalibrateCellTime (tcell4, GetCalorimeter(), absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
4337  }
4338 
4339  if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
4340  if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
4341  if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
4342  if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
4343 
4344  return ecell1+ecell2+ecell3+ecell4;
4345  }
4346  else // PHOS
4347  {
4348  Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
4349 
4350  Int_t relId1[] = { imod+1, 0, irow+1, icol };
4351  Int_t relId2[] = { imod+1, 0, irow-1, icol };
4352  Int_t relId3[] = { imod+1, 0, irow , icol+1 };
4353  Int_t relId4[] = { imod+1, 0, irow , icol-1 };
4354 
4355  GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
4356  GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
4357  GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
4358  GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
4359 
4360  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
4361 
4362  if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
4363  if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
4364  if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
4365  if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
4366 
4367  return ecell1+ecell2+ecell3+ecell4;
4368  }
4369 }
4370 
4371 //___________________________________________________________________________________________________________
4377 //___________________________________________________________________________________________________________
4378 void AliAnaCalorimeterQA::InvariantMassHistograms(Int_t iclus, Int_t nModule, const TObjArray* caloClusters,
4379  AliVCaloCells * cells)
4380 {
4381  AliDebug(1,"Start");
4382 
4383  //Get vertex for photon momentum calculation and event selection
4384  Double_t v[3] = {0,0,0}; //vertex ;
4385  //GetReader()->GetVertex(v);
4386 
4387  Int_t nModule2 = -1;
4388  Int_t nCaloClusters = caloClusters->GetEntriesFast();
4389 
4390  Float_t phi1 = fClusterMomentum.Phi();
4391  if(phi1 < 0) phi1 += TMath::TwoPi();
4392 
4393  Double_t tof1 = ((AliVCluster*) caloClusters->At(iclus))->GetTOF()*1.e9;
4394  if(tof1>400) tof1-=fConstantTimeShift;
4395 
4396  for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
4397  {
4398  AliVCluster* clus2 = (AliVCluster*) caloClusters->At(jclus);
4399 
4400  Float_t maxCellFraction = 0.;
4401  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2, maxCellFraction);
4402 
4403  Double_t tof2 = clus2->GetTOF()*1.e9;
4404  if(tof2>400) tof2-=fConstantTimeShift;
4405 
4406  Double_t diffTof = tof1-tof2;
4407 
4408  // Try to reduce background with a mild shower shape cut and no more
4409  // than 1 local maximum in cluster and remove low energy clusters
4410 
4411  if( !IsGoodCluster(absIdMax, clus2->GetM02(), clus2->GetNCells(), cells)
4412  || GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1
4413  || clus2->GetM02() > fInvMassMaxM02Cut
4414  || clus2->GetM02() < fInvMassMinM02Cut
4415  || clus2->E() < fInvMassMinECut
4416  || clus2->E() > fInvMassMaxECut
4417  || TMath::Abs(diffTof) > fInvMassMaxTimeDifference
4418  ) continue;
4419 
4420  // Get cluster kinematics
4421  clus2->GetMomentum(fClusterMomentum2,v);
4422 
4423  // Check only certain regions
4424  Bool_t in2 = kTRUE;
4425  if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(fClusterMomentum2.Eta(),fClusterMomentum2.Phi(),GetCalorimeter()) ;
4426  if(!in2) continue;
4427 
4428  Float_t pairPt = (fClusterMomentum+fClusterMomentum2).Pt();
4429 
4430  // Opening angle cut, avoid combination of DCal and EMCal clusters
4431  Double_t angle = fClusterMomentum.Angle(fClusterMomentum2.Vect());
4432 
4433  if ( fFillInvMassOpenAngle ) fhOpAngle->Fill(pairPt, angle, GetEventWeight()) ;
4434 
4435  if( angle > fInvMassMaxOpenAngle ) continue;
4436 
4437  // Get module of cluster
4438  nModule2 = GetModuleNumber(clus2);
4439 
4440  // Fill histograms
4441  Float_t mass = (fClusterMomentum+fClusterMomentum2).M ();
4442  Float_t asym = TMath::Abs( fClusterMomentum.E() - fClusterMomentum2.E() ) /
4443  ( fClusterMomentum.E() + fClusterMomentum2.E() );
4444 
4445  // All modules
4446  // Combine EMCal or PHOS clusters
4447 
4448  Float_t phi2 = fClusterMomentum2.Phi();
4449  if(phi2 < 0) phi2 += TMath::TwoPi();
4450 
4451  AliDebug(1,Form("Selected pair: pT %f, mass %f, asym %f, angle %f, diffTof %f, SM1 %d, SM2 %d\n",pairPt,mass,asym,angle,diffTof,nModule,nModule2));
4452 
4453  Bool_t inPi0Window = kFALSE;
4454  if(mass < 0.18 && mass > 0.1) inPi0Window = kTRUE ;
4455 
4456  if ( nModule < 12 && nModule2 < 12 )
4457  {
4458  fhIM ->Fill(pairPt, mass, GetEventWeight());
4459 
4460  if( fFillPi0PairDiffTime && inPi0Window )
4461  fhClusterPairDiffTimeEPi0Mass->Fill(pairPt, diffTof, GetEventWeight());
4462 
4463  if ( nModule == nModule2 )
4464  {
4465  fhIMSame->Fill(pairPt, mass, GetEventWeight());
4466 
4467  if( fFillPi0PairDiffTime && inPi0Window )
4468  fhClusterPairDiffTimeEPi0MassSame->Fill(pairPt, diffTof, GetEventWeight());
4469  }
4470  else
4471  {
4472  fhIMDiff->Fill(pairPt, mass, GetEventWeight());
4473  }
4474  }
4475  // Combine DCal clusters
4476  else if ( ( GetCalorimeter() == kEMCAL || GetCalorimeter() == kDCAL ) &&
4477  nModule > 11 && nModule2 > 11 && fNModules > 12 )
4478  {
4479  fhIMDCAL->Fill(pairPt, mass, GetEventWeight());
4480 
4481  if( fFillPi0PairDiffTime && inPi0Window )
4482  fhClusterPairDiffTimeEPi0MassDCal->Fill(pairPt, diffTof, GetEventWeight());
4483 
4484  if ( nModule == nModule2 )
4485  {
4486  fhIMDCALSame->Fill(pairPt, mass, GetEventWeight());
4487 
4488  if( fFillPi0PairDiffTime && inPi0Window )
4489  fhClusterPairDiffTimeEPi0MassDCalSame->Fill(pairPt, diffTof, GetEventWeight());
4490  }
4491  else
4492  {
4493  fhIMDCALDiff->Fill(pairPt, mass, GetEventWeight());
4494  }
4495  }
4496 
4497  if ( fFillInvMassOpenAngle ) fhIMvsOpAngle->Fill(mass, angle, GetEventWeight());
4498 
4499  // Single module
4500  if(nModule == nModule2 && nModule >= 0 && nModule < fNModules)
4501  fhIMMod[nModule]->Fill(pairPt, mass, GetEventWeight());
4502 
4503  // Asymetry histograms
4504  fhAsym->Fill(pairPt, asym, GetEventWeight());
4505  } // 2nd cluster loop
4506 
4507  //
4508  // Combine PHOS and DCal
4509  //
4510 
4511  if( ( GetCalorimeter() == kEMCAL || GetCalorimeter() == kDCAL ) &&
4512  fNModules > 12 && nModule > 11)
4513  {
4514  AliDebug(1,"Check DCal-PHOS pairs\n");
4515 
4516  Int_t sector1 = -1;
4517  Int_t sector2 = -1;
4518 
4519  if(phi1 >= 260*TMath::DegToRad() && phi1 < 280) sector1 = 0;
4520  if(phi1 >= 280*TMath::DegToRad() && phi1 < 300) sector1 = 1;
4521  if(phi1 >= 300*TMath::DegToRad() && phi1 < 320) sector1 = 2;
4522 
4523  for(Int_t jclus = 0 ; jclus < GetPHOSClusters()->GetEntriesFast() ; jclus++)
4524  {
4525  AliVCluster* clus2 = (AliVCluster*) GetPHOSClusters()->At(jclus);
4526 
4527  Float_t maxCellFraction = 0.;
4528  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2, maxCellFraction);
4529 
4530  // Try to reduce background, remove low energy clusters
4531  if( !IsGoodCluster(absIdMax, clus2->GetM02(), clus2->GetNCells(), cells)
4532  || clus2->E() < fInvMassMinECut
4533  || clus2->E() > fInvMassMaxECut
4534  ) continue;
4535 
4536  // Get cluster kinematics
4537  clus2->GetMomentum(fClusterMomentum2,v);
4538 
4539  // Fill histograms
4540 
4541  Float_t mass = (fClusterMomentum+fClusterMomentum2).M ();
4542  Float_t pairPt = (fClusterMomentum+fClusterMomentum2).Pt();
4543  //Float_t asym = TMath::Abs( fClusterMomentum.E() - fClusterMomentum2.E() ) /
4544  //( fClusterMomentum.E() + fClusterMomentum2.E() );
4545 
4546  fhIMDCALPHOS->Fill(pairPt, mass, GetEventWeight());
4547 
4548  Float_t phiPHOS = fClusterMomentum2.Phi();
4549  if(phiPHOS < 0) phiPHOS += TMath::TwoPi();
4550 
4551  if(phiPHOS >= 260*TMath::DegToRad() && phiPHOS < 280) sector2 = 0;
4552  if(phiPHOS >= 280*TMath::DegToRad() && phiPHOS < 300) sector2 = 1;
4553  if(phiPHOS >= 300*TMath::DegToRad() && phiPHOS < 320) sector2 = 2;
4554 
4555  if(sector1 == sector2) fhIMDCALPHOSSame->Fill(pairPt, mass, GetEventWeight());
4556 
4557  } // PHOS cluster loop
4558  } // DCal-PHOS combination
4559 
4560 
4561  // Select EMCal clusters in DCal eta acceptance
4562  // Cross check combination of DCal-PHOS pairs
4563  if( fFillInvMassInEMCALWithPHOSDCalAcc && TMath::Abs(fClusterMomentum.Eta() > 0.22))
4564  {
4565  AliDebug(1,"Check EMCAL(DCal)-EMCAL(PHOS) pairs\n");
4566 
4567  Int_t sector1 = -1;
4568  Int_t sector2 = -1;
4569 
4570  if(phi1 >= 80*TMath::DegToRad() && phi1 < 100) sector1 = 0;
4571  if(phi1 >= 100*TMath::DegToRad() && phi1 < 120) sector1 = 1;
4572  if(phi1 >= 120*TMath::DegToRad() && phi1 < 140) sector1 = 2;
4573  if(phi1 >= 140*TMath::DegToRad() && phi1 < 160) sector1 = 3;
4574  if(phi1 >= 160*TMath::DegToRad() && phi1 < 180) sector1 = 4;
4575  if(phi1 >= 180*TMath::DegToRad() && phi1 < 190) sector1 = 5;
4576 
4577  for(Int_t jclus = 0 ; jclus < caloClusters->GetEntriesFast() ; jclus++)
4578  {
4579  AliVCluster* clus2 = (AliVCluster*) caloClusters->At(jclus);
4580 
4581  Float_t maxCellFraction = 0.;
4582  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2, maxCellFraction);
4583 
4584  Double_t tof2 = clus2->GetTOF()*1.e9;
4585  if(tof2>400) tof2-=fConstantTimeShift;
4586 
4587  Double_t diffTof = TMath::Abs(tof1-tof2);
4588 
4589  // Try to reduce background with a mild shower shape cut and no more
4590  // than 1 local maximum in cluster and remove low energy clusters
4591  if( !IsGoodCluster(absIdMax, clus2->GetM02(), clus2->GetNCells(), cells)
4592  || GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1
4593  || clus2->GetM02() > fInvMassMaxM02Cut
4594  || clus2->GetM02() < fInvMassMinM02Cut
4595  || clus2->E() < fInvMassMinECut
4596  || clus2->E() > fInvMassMaxECut
4597  || TMath::Abs(diffTof) > fInvMassMaxTimeDifference
4598  ) continue;
4599 
4600  // Get cluster kinematics
4601  clus2->GetMomentum(fClusterMomentum2,v);
4602 
4603  // Select EMCal clusters in PHOS acceptance
4604  if(TMath::Abs(fClusterMomentum2.Eta() > 0.13)) continue ;
4605 
4606  // Fill histograms
4607 
4608  Float_t mass = (fClusterMomentum+fClusterMomentum2).M ();
4609  Float_t pairPt = (fClusterMomentum+fClusterMomentum2).Pt();
4610  //Float_t asym = TMath::Abs( fClusterMomentum.E() - fClusterMomentum2.E() ) /
4611  //( fClusterMomentum.E() + fClusterMomentum2.E() );
4612 
4613  fhIMEMCALPHOS->Fill(pairPt, mass, GetEventWeight());
4614 
4615  Float_t phiPHOS = fClusterMomentum2.Phi();
4616  if(phiPHOS < 0) phiPHOS += TMath::TwoPi();
4617 
4618  if(phiPHOS >= 80*TMath::DegToRad() && phiPHOS < 100) sector2 = 0;
4619  if(phiPHOS >= 100*TMath::DegToRad() && phiPHOS < 120) sector2 = 1;
4620  if(phiPHOS >= 120*TMath::DegToRad() && phiPHOS < 140) sector2 = 2;
4621  if(phiPHOS >= 140*TMath::DegToRad() && phiPHOS < 160) sector2 = 3;
4622  if(phiPHOS >= 160*TMath::DegToRad() && phiPHOS < 180) sector2 = 4;
4623  if(phiPHOS >= 180*TMath::DegToRad() && phiPHOS < 190) sector2 = 5;
4624 
4625  if(sector1 == sector2) fhIMEMCALPHOSSame->Fill(pairPt, mass, GetEventWeight());
4626 
4627  } // PHOS cluster loop
4628  } // EMCal(DCal)-EMCAL(PHOS) combination
4629 
4630  AliDebug(1,"End");
4631 }
4632 
4633 //______________________________
4634 // Check if the calorimeter setting is ok, if not abort.
4635 //______________________________
4636 void AliAnaCalorimeterQA::Init()
4637 {
4638  if(GetCalorimeter() != kPHOS && GetCalorimeter() !=kEMCAL)
4639  AliFatal(Form("Wrong calorimeter name <%s>", GetCalorimeterString().Data()));
4640 
4641  //if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
4642  // AliFatal("Analysis of reconstructed data, MC reader not aplicable");
4643 }
4644 
4645 //________________________________________
4647 //________________________________________
4648 void AliAnaCalorimeterQA::InitParameters()
4649 {
4650  AddToHistogramsName("AnaCaloQA_");
4651 
4652  fNModules = 22; // set maximum to maximum number of EMCAL modules
4653  fNRCU = 2; // set maximum number of RCU in EMCAL per SM
4654 
4655  fTimeCutMin = -9999999;
4656  fTimeCutMax = 9999999;
4657 
4658  fEMCALCellAmpMin = 0.2; // 200 MeV
4659  fPHOSCellAmpMin = 0.2; // 200 MeV
4660  fCellAmpMin = 0.2; // 200 MeV
4661 
4662  fEMCALClusterNCellMin = 2; // at least 2
4663  fPHOSClusterNCellMin = 3; // at least 3
4664 
4665  fInvMassMinECut = 0.5; // 500 MeV
4666  fInvMassMaxECut = 10;
4667 
4668  fInvMassMinM02Cut = 0.1;
4669  fInvMassMaxM02Cut = 0.4;
4670 
4671  fInvMassMaxTimeDifference = 200; // ns, quite open
4672 
4673  fInvMassMaxOpenAngle = 100*TMath::DegToRad(); // 100 degrees
4674 
4675  // Exotic studies
4676  fExoNECrossCuts = 10 ;
4677  fExoNDTimeCuts = 4 ;
4678 
4679  fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
4680  fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
4681  fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
4682 
4683  fNEBinCuts = 14;
4684  fEBinCuts[0] = 0.; fEBinCuts[1] = 0.3; fEBinCuts[2] = 0.5;
4685  fEBinCuts[3] = 1.; fEBinCuts[4] = 2. ; fEBinCuts[5] = 3. ;
4686  fEBinCuts[6] = 4.; fEBinCuts[7] = 5. ; fEBinCuts[8] = 7. ;
4687  fEBinCuts[9] = 9.; fEBinCuts[10]= 12.; fEBinCuts[11]= 15.;
4688  fEBinCuts[12]= 20.; fEBinCuts[13]= 50.; fEBinCuts[14]= 100.;
4689  for(Int_t i = fNEBinCuts; i < 15; i++) fEBinCuts[i] = 1000.;
4690 
4691 
4692 }
4693 
4694 //_____________________________________________________________________________
4696 //_____________________________________________________________________________
4697 Bool_t AliAnaCalorimeterQA::IsGoodCluster(Int_t absIdMax, Float_t m02,
4698  Int_t nCellsPerCluster, AliVCaloCells* cells)
4699 {
4700  //if(!fStudyBadClusters) return kTRUE;
4701 
4702  if(GetCalorimeter() == kEMCAL)
4703  {
4704 
4705  if( m02 < 0.05 || nCellsPerCluster < fEMCALClusterNCellMin ) return kFALSE ; // mild shower shape cut for exotics
4706 
4707  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
4708  {
4709  return !(GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()));
4710  }
4711  else
4712  {
4713  return kTRUE;
4714  }
4715  }
4716  else // PHOS
4717  {
4718  if( m02 < 0.05 || nCellsPerCluster < fPHOSClusterNCellMin ) return kFALSE ; // mild shower shape cut for exotics
4719 
4720  Float_t ampMax = cells->GetCellAmplitude(absIdMax);
4721  GetCaloUtils()->RecalibrateCellAmplitude(ampMax, GetCalorimeter(), absIdMax);
4722