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