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