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