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