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