]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx
Allow fill QA time histograms for Cells in AODs, add a switch to fill them
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaCalorimeterQA.cxx
CommitLineData
9725fd2a 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 **************************************************************************/
9725fd2a 15
a6f26052 16//_________________________________________________________________________
17// Class to check results from simulations or reconstructed real data.
18// Fill few histograms and do some checking plots
19//
20//-- Author: Gustavo Conesa (INFN-LNF)
21//_________________________________________________________________________
9725fd2a 22
23
a6f26052 24// --- ROOT system ---
d55bb5e1 25#include <TObjArray.h>
26#include <TParticle.h>
27#include <TDatabasePDG.h>
28#include <TH3F.h>
0c1383b5 29#include <TObjString.h>
9725fd2a 30
a6f26052 31//---- AliRoot system ----
9725fd2a 32#include "AliAnaCalorimeterQA.h"
33#include "AliCaloTrackReader.h"
34#include "AliStack.h"
c8fe2783 35#include "AliVCaloCells.h"
ff45398a 36#include "AliFiducialCut.h"
c8fe2783 37#include "AliVCluster.h"
d55bb5e1 38#include "AliVTrack.h"
c8fe2783 39#include "AliVEvent.h"
902aa95c 40#include "AliVEventHandler.h"
902aa95c 41#include "AliAODMCParticle.h"
42#include "AliMCAnalysisUtils.h"
9725fd2a 43
c5693f62 44// --- Detectors ---
45#include "AliPHOSGeoUtils.h"
46#include "AliEMCALGeometry.h"
47
9725fd2a 48ClassImp(AliAnaCalorimeterQA)
c8fe2783 49
649b825d 50//________________________________________
c8fe2783 51AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
745913ae 52AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""),
649b825d 53
54//Switches
e6fec6f5 55fFillAllCellTimeHisto(kTRUE),
9e9f04cb 56fFillAllPosHisto(kFALSE), fFillAllPosHisto2(kTRUE),
57fFillAllTH12(kFALSE), fFillAllTH3(kTRUE),
35c71d5c 58fFillAllTMHisto(kTRUE), fFillAllPi0Histo(kTRUE),
649b825d 59fCorrelate(kTRUE), fStudyBadClusters(kFALSE),
60fStudyClustersAsymmetry(kFALSE), fStudyWeight(kFALSE),
61
62//Parameters and cuts
35c71d5c 63fNModules(12), fNRCU(2),
64fNMaxCols(48), fNMaxRows(24),
f15c25da 65fTimeCutMin(-10000), fTimeCutMax(10000),
9e9f04cb 66fEMCALCellAmpMin(0), fPHOSCellAmpMin(0),
649b825d 67
68//Histograms
9e9f04cb 69fhE(0), fhPt(0),
70fhPhi(0), fhEta(0), fhEtaPhiE(0),
71fhECharged(0), fhPtCharged(0),
72fhPhiCharged(0), fhEtaCharged(0), fhEtaPhiECharged(0),
521636d2 73
74//Invariant mass
9e9f04cb 75fhIM(0 ), fhAsym(0),
a82b4462 76
77fhNCellsPerCluster(0), fhNCellsPerClusterNoCut(0), fhNClusters(0),
3129a79e 78
e1e62b89 79//Timing
9e9f04cb 80fhClusterTimeEnergy(0), fhCellTimeSpreadRespectToCellMax(0),
9e9f04cb 81fhCellIdCellLargeTimeSpread(0), fhClusterPairDiffTimeE(0),
9e9f04cb 82fhClusterMaxCellCloseCellRatio(0), fhClusterMaxCellCloseCellDiff(0),
83fhClusterMaxCellDiff(0), fhClusterMaxCellDiffNoCut(0),
a82b4462 84fhClusterMaxCellDiffAverageTime(0), fhClusterMaxCellDiffWeightedTime(0),
1a72f6c5 85fhClusterMaxCellECross(0),
649b825d 86fhLambda0(0), fhLambda1(0), fhDispersion(0),
715fd81f 87
649b825d 88//bad clusters
89fhBadClusterEnergy(0), fhBadClusterTimeEnergy(0),
90fhBadClusterPairDiffTimeE(0), fhBadCellTimeSpreadRespectToCellMax(0),
9e9f04cb 91fhBadClusterMaxCellCloseCellRatio(0), fhBadClusterMaxCellCloseCellDiff(0), fhBadClusterMaxCellDiff(0),
a82b4462 92fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
1a72f6c5 93fhBadClusterMaxCellECross(0),
1a83b960 94fhBadClusterDeltaIEtaDeltaIPhiE0(0), fhBadClusterDeltaIEtaDeltaIPhiE2(0),
95fhBadClusterDeltaIEtaDeltaIPhiE6(0), fhBadClusterDeltaIA(0),
9e9f04cb 96
521636d2 97//Position
9e9f04cb 98fhRNCells(0), fhXNCells(0),
99fhYNCells(0), fhZNCells(0),
100fhRE(0), fhXE(0),
101fhYE(0), fhZE(0),
521636d2 102fhXYZ(0),
9e9f04cb 103fhRCellE(0), fhXCellE(0),
104fhYCellE(0), fhZCellE(0),
521636d2 105fhXYZCell(0),
9e9f04cb 106fhDeltaCellClusterRNCells(0), fhDeltaCellClusterXNCells(0),
107fhDeltaCellClusterYNCells(0), fhDeltaCellClusterZNCells(0),
108fhDeltaCellClusterRE(0), fhDeltaCellClusterXE(0),
109fhDeltaCellClusterYE(0), fhDeltaCellClusterZE(0),
649b825d 110
521636d2 111// Cells
9e9f04cb 112fhNCells(0), fhAmplitude(0),
113fhAmpId(0), fhEtaPhiAmp(0),
1a72f6c5 114fhTime(0), fhTimeVz(0),
115fhTimeId(0), fhTimeAmp(0),
116fhCellECross(0),
9e9f04cb 117fhCaloCorrNClusters(0), fhCaloCorrEClusters(0),
118fhCaloCorrNCells(0), fhCaloCorrECells(0),
119fhCaloV0SCorrNClusters(0), fhCaloV0SCorrEClusters(0),
120fhCaloV0SCorrNCells(0), fhCaloV0SCorrECells(0),
121fhCaloV0MCorrNClusters(0), fhCaloV0MCorrEClusters(0),
122fhCaloV0MCorrNCells(0), fhCaloV0MCorrECells(0),
123fhCaloTrackMCorrNClusters(0), fhCaloTrackMCorrEClusters(0),
124fhCaloTrackMCorrNCells(0), fhCaloTrackMCorrECells(0),
649b825d 125
521636d2 126//Super-Module dependent histgrams
649b825d 127fhEMod(0), fhAmpMod(0), fhTimeMod(0),
128fhNClustersMod(0), fhNCellsMod(0),
129fhNCellsPerClusterMod(0), fhNCellsPerClusterModNoCut(0),
130
131fhGridCells(0), fhGridCellsE(0), fhGridCellsTime(0),
132fhTimeAmpPerRCU(0), fhIMMod(0),
133
134// Weight studies
135fhECellClusterRatio(0), fhECellClusterLogRatio(0),
136fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
715fd81f 137
138// MC and reco
649b825d 139fhRecoMCE(), fhRecoMCPhi(), fhRecoMCEta(),
140fhRecoMCDeltaE(), fhRecoMCRatioE(),
141fhRecoMCDeltaPhi(), fhRecoMCDeltaEta(),
715fd81f 142
521636d2 143// MC only
649b825d 144fhGenMCE(), fhGenMCEtaPhi(),
145fhGenMCAccE(), fhGenMCAccEtaPhi(),
35c71d5c 146
147//matched MC
649b825d 148fhEMVxyz(0), fhEMR(0),
149fhHaVxyz(0), fhHaR(0),
d55bb5e1 150fh1EOverP(0), fh2dR(0),
649b825d 151fh2EledEdx(0), fh2MatchdEdx(0),
d55bb5e1 152fhMCEle1EOverP(0), fhMCEle1dR(0), fhMCEle2MatchdEdx(0),
153fhMCChHad1EOverP(0), fhMCChHad1dR(0), fhMCChHad2MatchdEdx(0),
154fhMCNeutral1EOverP(0), fhMCNeutral1dR(0), fhMCNeutral2MatchdEdx(0), fh1EOverPR02(0),
155fhMCEle1EOverPR02(0), fhMCChHad1EOverPR02(0), fhMCNeutral1EOverPR02(0)
9725fd2a 156{
a6f26052 157 //Default Ctor
afabc52f 158
649b825d 159 //Weight studies
1a72f6c5 160 for(Int_t i =0; i < 14; i++){
649b825d 161 fhLambda0ForW0[i] = 0;
1a72f6c5 162 //fhLambda1ForW0[i] = 0;
649b825d 163
164 for(Int_t j = 0; j < 5; j++){
165 fhLambda0ForW0MC[i][j] = 0;
1a72f6c5 166 //fhLambda1ForW0MC[i][j] = 0;
649b825d 167 }
168
169 }
c8fe2783 170
649b825d 171 //Cluster size
172 fhDeltaIEtaDeltaIPhiE0[0] = 0 ; fhDeltaIEtaDeltaIPhiE2[0] = 0; fhDeltaIEtaDeltaIPhiE6[0] = 0;
173 fhDeltaIEtaDeltaIPhiE0[1] = 0 ; fhDeltaIEtaDeltaIPhiE2[1] = 0; fhDeltaIEtaDeltaIPhiE6[1] = 0;
174 fhDeltaIA[0] = 0 ; fhDeltaIAL0[0] = 0; fhDeltaIAL1[0] = 0;
175 fhDeltaIA[1] = 0 ; fhDeltaIAL0[1] = 0; fhDeltaIAL1[1] = 0;
176 fhDeltaIANCells[0] = 0 ; fhDeltaIANCells[1] = 0;
177 fhDeltaIAMC[0] = 0 ; fhDeltaIAMC[1] = 0;
178 fhDeltaIAMC[2] = 0 ; fhDeltaIAMC[3] = 0;
2302a644 179
649b825d 180 // MC
c8fe2783 181
649b825d 182 for(Int_t i = 0; i < 6; i++){
183
184 fhRecoMCE[i][0] = 0; fhRecoMCE[i][1] = 0;
185 fhRecoMCPhi[i][0] = 0; fhRecoMCPhi[i][1] = 0;
186 fhRecoMCEta[i][0] = 0; fhRecoMCEta[i][1] = 0;
187 fhRecoMCDeltaE[i][0] = 0; fhRecoMCDeltaE[i][1] = 0;
188 fhRecoMCRatioE[i][0] = 0; fhRecoMCRatioE[i][1] = 0;
189 fhRecoMCDeltaPhi[i][0] = 0; fhRecoMCDeltaPhi[i][1] = 0;
190 fhRecoMCDeltaEta[i][0] = 0; fhRecoMCDeltaEta[i][1] = 0;
191
192 }
2302a644 193
649b825d 194 //Initialize parameters
195 InitParameters();
196}
3f5990d6 197
1a83b960 198//____________________________________________________________________________________________________________________________
c5693f62 199void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
1a83b960 200 const Int_t absIdMax, const Double_t maxCellFraction, const Float_t eCrossFrac,
a82b4462 201 const Double_t tmax, Double_t timeAverages[2]
649b825d 202 )
203{
204 //Bad cluster histograms
b0114dba 205
206 // printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
207 // GetReader()->GetEventNumber(), fCalorimeter.Data(),
208 // clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
1a72f6c5 209
649b825d 210 fhBadClusterEnergy ->Fill(clus->E());
211 Double_t tof = clus->GetTOF()*1.e9;
1a83b960 212 fhBadClusterTimeEnergy ->Fill(clus->E(),tof);
213 fhBadClusterMaxCellDiff ->Fill(clus->E(),maxCellFraction);
214 fhBadClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
215
216 if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
1a72f6c5 217
a82b4462 218 //Clusters in event time differencem bad minus good
2302a644 219
649b825d 220 for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
a6f26052 221
649b825d 222 AliVCluster* clus2 = (AliVCluster*)caloClusters->At(iclus2);
a6f26052 223
649b825d 224 if(clus->GetID()==clus2->GetID()) continue;
715fd81f 225
a82b4462 226 Float_t maxCellFraction2 = 0.;
227 Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
228 if(IsGoodCluster(absIdMax2,cells)){
229 Double_t tof2 = clus2->GetTOF()*1.e9;
649b825d 230 fhBadClusterPairDiffTimeE ->Fill(clus->E(), (tof-tof2));
649b825d 231 }
a82b4462 232
649b825d 233 } // loop
924e319f 234
649b825d 235 // Max cell compared to other cells in cluster
e6fec6f5 236 if(fFillAllCellTimeHisto)
2747966a 237 {
649b825d 238 fhBadClusterMaxCellDiffAverageTime ->Fill(clus->E(),tmax-timeAverages[0]);
649b825d 239 fhBadClusterMaxCellDiffWeightedTime ->Fill(clus->E(),tmax-timeAverages[1]);
649b825d 240 }
715fd81f 241
2747966a 242 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
243 {
649b825d 244 Int_t absId = clus->GetCellsAbsId()[ipos];
2747966a 245 if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
246 {
649b825d 247 Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
248
249 fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
250 fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
251
e6fec6f5 252 if(fFillAllCellTimeHisto)
dbba06ca 253 {
649b825d 254 Double_t time = cells->GetCellTime(absId);
dbba06ca 255 GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
256
649b825d 257 Float_t diff = (tmax-time*1e9);
258 fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
259
e6fec6f5 260 }
649b825d 261 }// Not max
262 }//loop
715fd81f 263
649b825d 264}
715fd81f 265
dbf54f1e 266//______________________________________________________________________
267void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus,
268 AliVCaloCells* cells,
269 Double_t timeAverages[2])
649b825d 270{
271 // Calculate time averages and weights
a95eac90 272
649b825d 273 // First recalculate energy in case non linearity was applied
274 Float_t energy = 0;
275 Float_t ampMax = 0, amp = 0;
276 Int_t absIdMax =-1;
2747966a 277 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
278 {
649b825d 279 Int_t id = clus->GetCellsAbsId()[ipos];
e1e62b89 280
649b825d 281 //Recalibrate cell energy if needed
282 amp = cells->GetCellAmplitude(id);
dbba06ca 283 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
4c8f7c2e 284
649b825d 285 energy += amp;
e1e62b89 286
2747966a 287 if(amp> ampMax)
288 {
649b825d 289 ampMax = amp;
290 absIdMax = id;
291 }
4c8f7c2e 292
649b825d 293 } // energy loop
294
295 // Calculate average time of cells in cluster and weighted average
dbf54f1e 296 Double_t aTime = 0;
297 Double_t wTime = 0;
298 Float_t wTot = 0;
299 Double_t time = 0;
300 Int_t id =-1;
301 Double_t w = 0;
302 Int_t ncells = clus->GetNCells();
2747966a 303 for (Int_t ipos = 0; ipos < ncells; ipos++)
304 {
dbf54f1e 305 id = clus ->GetCellsAbsId()[ipos];
306 amp = cells->GetCellAmplitude(id);
307 time = cells->GetCellTime(id);
649b825d 308
309 //Recalibrate energy and time
dbba06ca 310 GetCaloUtils()->RecalibrateCellAmplitude(amp , fCalorimeter, id);
311 GetCaloUtils()->RecalibrateCellTime (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
312
dbf54f1e 313 w = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
649b825d 314 aTime += time*1e9;
315 wTime += time*1e9 * w;
316 wTot += w;
35c71d5c 317
649b825d 318 }
f3138ecf 319
dbf54f1e 320 if(ncells > 0) aTime /= ncells;
321 else aTime = 0;
649b825d 322
dbf54f1e 323 if(wTot > 0) wTime /= wTot;
f3138ecf 324 else wTime = 0;
325
dbf54f1e 326 timeAverages[0] = aTime;
327 timeAverages[1] = wTime;
39de6caa 328
649b825d 329}
330
331//____________________________________________________________
332void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
333{
334 // Plot histograms related to cells only
9e9f04cb 335
649b825d 336 Int_t ncells = cells->GetNumberOfCells();
39de6caa 337
649b825d 338 if(GetDebug() > 0)
339 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), ncells );
340
341 //Init arrays and used variables
342 Int_t *nCellsInModule = new Int_t[fNModules];
343 for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
a82b4462 344
649b825d 345 Int_t icol = -1;
346 Int_t irow = -1;
347 Int_t iRCU = -1;
348 Float_t amp = 0.;
349 Double_t time = 0.;
350 Int_t id = -1;
351 Float_t recalF = 1.;
a82b4462 352 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
649b825d 353
354 for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++) {
355 if(GetDebug() > 2)
356 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell));
357 Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),fCalorimeter, icol, irow, iRCU);
358 if(GetDebug() > 2)
359 printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
360
361 if(nModule < fNModules) {
362
363 //Check if the cell is a bad channel
364 if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn()){
ae2c2bc4 365 if(fCalorimeter=="EMCAL")
366 {
649b825d 367 if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
368 }
ae2c2bc4 369 else
370 {
371 if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
649b825d 372 }
373 } // use bad channel map
374
375 amp = cells->GetAmplitude(iCell)*recalF;
376 time = cells->GetTime(iCell);
377 id = cells->GetCellNumber(iCell);
378
379 // Amplitude recalibration if set
dbba06ca 380 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 381
382 // Time recalibration if set
dbba06ca 383 GetCaloUtils()->RecalibrateCellTime (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
649b825d 384
385 //Transform time to ns
386 time *= 1.0e9;
f15c25da 387
e6fec6f5 388 if(time < fTimeCutMin || time > fTimeCutMax)
a87e069d 389 {
649b825d 390 if(GetDebug() > 0 )printf("AliAnaCalorimeterQA - Remove cell with Time %f\n",time);
391 continue;
e6fec6f5 392 }
06f1b12a 393
394 // Remove exotic cells, defined only for EMCAL
395 if(fCalorimeter=="EMCAL" &&
396 GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
397
649b825d 398
399 fhAmplitude->Fill(amp);
400 fhAmpId ->Fill(amp,id);
401 fhAmpMod ->Fill(amp,nModule);
402
403 if ((fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ||
2747966a 404 (fCalorimeter=="PHOS" && amp > fPHOSCellAmpMin ) )
405 {
406
407 //E cross for exotic cells
408 if(amp > 0.01) fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
e3300002 409
649b825d 410 nCellsInModule[nModule]++ ;
06f1b12a 411
649b825d 412 Int_t icols = icol;
413 Int_t irows = irow;
57d8227a 414
415 if(fCalorimeter=="EMCAL")
416 {
649b825d 417 icols = (nModule % 2) ? icol + fNMaxCols : icol;
57d8227a 418 if(nModule < 10 )
419 irows = irow + fNMaxRows * Int_t(nModule / 2);
420 else // 1/3 SM
421 irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
649b825d 422 }
57d8227a 423 else
424 {
06f1b12a 425 irows = irow + fNMaxRows * nModule;
649b825d 426 }
06f1b12a 427
649b825d 428 fhGridCells ->Fill(icols,irows);
429 fhGridCellsE->Fill(icols,irows,amp);
430
e6fec6f5 431 if(fFillAllCellTimeHisto)
a87e069d 432 {
649b825d 433 //printf("%s: time %g\n",fCalorimeter.Data(), time);
1a72f6c5 434
435 Double_t v[3] = {0,0,0}; //vertex ;
436 GetReader()->GetVertex(v);
437 if(amp > 0.5) fhTimeVz ->Fill(TMath::Abs(v[2]),time);
438
649b825d 439 fhTime ->Fill(time);
440 fhTimeId ->Fill(time,id);
441 fhTimeAmp ->Fill(amp,time);
442 fhGridCellsTime->Fill(icols,irows,time);
443 fhTimeMod ->Fill(time,nModule);
444 fhTimeAmpPerRCU [nModule*fNRCU+iRCU]->Fill(amp, time);
445
446 }
447 }
448
449 //Get Eta-Phi position of Cell
450 if(fFillAllPosHisto)
451 {
452 if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
453 Float_t celleta = 0.;
454 Float_t cellphi = 0.;
455 GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
456
457 fhEtaPhiAmp->Fill(celleta,cellphi,amp);
458 Double_t cellpos[] = {0, 0, 0};
459 GetEMCALGeometry()->GetGlobal(id, cellpos);
460 fhXCellE->Fill(cellpos[0],amp) ;
461 fhYCellE->Fill(cellpos[1],amp) ;
462 fhZCellE->Fill(cellpos[2],amp) ;
463 Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
464 fhRCellE->Fill(rcell,amp) ;
465 fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2]) ;
466 }//EMCAL Cells
467 else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
468 TVector3 xyz;
469 Int_t relId[4], module;
470 Float_t xCell, zCell;
471
472 GetPHOSGeometry()->AbsToRelNumbering(id,relId);
473 module = relId[0];
474 GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
475 GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
476 Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
477 fhXCellE ->Fill(xyz.X(),amp) ;
478 fhYCellE ->Fill(xyz.Y(),amp) ;
479 fhZCellE ->Fill(xyz.Z(),amp) ;
480 fhRCellE ->Fill(rcell ,amp) ;
481 fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z()) ;
482 }//PHOS cells
483 }//fill cell position histograms
484
485 if (fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ncells ++ ;
486 else if(fCalorimeter=="PHOS" && amp > fPHOSCellAmpMin) ncells ++ ;
487 //else
488 // printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - no %s CELLS passed the analysis cut\n",fCalorimeter.Data());
489 }//nmodules
490 }//cell loop
491
492 if(ncells > 0 )fhNCells->Fill(ncells) ; //fill the cells after the cut
493
494 //Number of cells per module
495 for(Int_t imod = 0; imod < fNModules; imod++ ) {
496
497 if(GetDebug() > 1)
498 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]);
499
500 fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
501
502 }
503
504 delete [] nCellsInModule;
505
506}
507
508//__________________________________________________________________________
509void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
510{
511 // Fill histograms releated to cell position
512
513
514 Int_t nCaloCellsPerCluster = clus->GetNCells();
515 UShort_t * indexList = clus->GetCellsAbsId();
516 Float_t pos[3];
517 clus->GetPosition(pos);
518 Float_t clEnergy = clus->E();
519
520 //Loop on cluster cells
521 for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
522
523 // printf("Index %d\n",ipos);
524 Int_t absId = indexList[ipos];
525
526 //Get position of cell compare to cluster
527
528 if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
529
530 Double_t cellpos[] = {0, 0, 0};
531 GetEMCALGeometry()->GetGlobal(absId, cellpos);
532
533 fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ;
534 fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ;
535 fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
536
537 fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy) ;
538 fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy) ;
539 fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy) ;
540
541 Float_t r = TMath::Sqrt(pos[0] *pos[0] + pos[1] * pos[1] );
542 Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
543
544 fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
545 fhDeltaCellClusterRE ->Fill(r-rcell, clEnergy) ;
546
547 }//EMCAL and its matrices are available
548 else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
549 TVector3 xyz;
550 Int_t relId[4], module;
551 Float_t xCell, zCell;
552
553 GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
554 module = relId[0];
555 GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
556 GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
557
558 fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ;
559 fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ;
560 fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
561
562 fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy) ;
563 fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy) ;
564 fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy) ;
565
566 Float_t r = TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1] );
567 Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
568
569 fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
570 fhDeltaCellClusterRE ->Fill(r-rcell, clEnergy) ;
571
572 }//PHOS and its matrices are available
573 }// cluster cell loop
574}
575
576//___________________________________________________________________________________________
1a83b960 577void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, const Int_t absIdMax,
578 const Bool_t goodCluster)
649b825d 579{
580 // Study the shape of the cluster in cell units terms
581
582 //No use to study clusters with less than 4 cells
583 if(clus->GetNCells() <=3 ) return;
584
585 Int_t dIeta = 0;
586 Int_t dIphi = 0;
587
588 Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
589 Int_t smMax = GetModuleNumberCellIndexes(absIdMax,fCalorimeter, ietaMax, iphiMax, rcuMax);
590
591 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
592
593 Int_t absId = clus->GetCellsAbsId()[ipos];
594
595 Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
596 Int_t sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ieta, iphi, rcu);
597
598 if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
599
600 if(smMax==sm){
601 if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
602 }
603 else {
604 Int_t ietaShift = ieta;
605 Int_t ietaMaxShift = ietaMax;
606 if (ieta > ietaMax) ietaMaxShift+=48;
607 else ietaShift +=48;
608 if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
609 }
610
649b825d 611 }// fill cell-cluster histogram loop
612
649b825d 613
1a83b960 614 Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
615
616 if(goodCluster)
617 {
649b825d 618
1a83b960 619 // Was cluster matched?
620 Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
649b825d 621
1a83b960 622 if (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
623 else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
624 else fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
625
626 fhDeltaIA[matched]->Fill(clus->E(),dIA);
627
628 if(clus->E() > 0.5){
629
630 fhDeltaIAL0[matched] ->Fill(clus->GetM02(),dIA);
631 fhDeltaIAL1[matched] ->Fill(clus->GetM20(),dIA);
632 fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
633
649b825d 634 }
635
1a83b960 636 // Origin of clusters
637 Int_t nLabel = clus->GetNLabels();
638 Int_t* labels = clus->GetLabels();
639 if(IsDataMC()){
640 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),0);
641 if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
642 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) &&
643 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
644 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
645 fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
646 }
647 else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
648 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
649 fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
650 }
651 else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
652 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
653 fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
654 }
655 else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
656 fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
657 }
658
659 } // MC
660 } // good cluster
661 else
662 {
663 if (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
664 else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
665 else fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
666
667 fhBadClusterDeltaIA->Fill(clus->E(),dIA);
668
669 }
649b825d 670}
671
1a83b960 672//_________________________________________________________________________________________________________________________
c5693f62 673void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus,const TObjArray *caloClusters, AliVCaloCells * cells,
1a83b960 674 const Int_t absIdMax, const Double_t maxCellFraction, const Float_t eCrossFrac,
a82b4462 675 const Double_t tmax, Double_t timeAverages[2])
649b825d 676{
677 //Fill CaloCluster related histograms
678
1a83b960 679 Double_t tof = clus->GetTOF()*1.e9;
649b825d 680
1a83b960 681 fhLambda0 ->Fill(clus->E(),clus->GetM02());
682 fhLambda1 ->Fill(clus->E(),clus->GetM20());
683 fhDispersion ->Fill(clus->E(),clus->GetDispersion());
649b825d 684
a82b4462 685 fhClusterMaxCellDiff ->Fill(clus->E(),maxCellFraction);
686 fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
687 fhClusterTimeEnergy ->Fill(clus->E(),tof);
649b825d 688
1a83b960 689 if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
690
649b825d 691 //Clusters in event time difference
692 for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
693
694 AliVCluster* clus2 = (AliVCluster*) caloClusters->At(iclus2);
695
696 if(clus->GetID()==clus2->GetID()) continue;
697
a87e069d 698 if(clus->GetM02() > 0.01 && clus2->GetM02() > 0.01)
699 {
a82b4462 700 Double_t tof2 = clus2->GetTOF()*1.e9;
649b825d 701 fhClusterPairDiffTimeE ->Fill(clus->E(), tof-tof2);
702 }
703 }
704
1a83b960 705 Int_t nModule = GetModuleNumber(clus);
706 Int_t nCaloCellsPerCluster = clus->GetNCells();
707
649b825d 708 if(nCaloCellsPerCluster > 1){
709
710 // check time of cells respect to max energy cell
711
e6fec6f5 712 if(fFillAllCellTimeHisto)
a87e069d 713 {
649b825d 714 fhClusterMaxCellDiffAverageTime ->Fill(clus->E(),tmax-timeAverages[0]);
649b825d 715 fhClusterMaxCellDiffWeightedTime ->Fill(clus->E(),tmax-timeAverages[1]);
649b825d 716 }
717
dbba06ca 718 for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
719 {
649b825d 720 Int_t absId = clus->GetCellsAbsId()[ipos];
2747966a 721 if(absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01) continue;
649b825d 722
723 Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
724 fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
725 fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
726
e6fec6f5 727 if(fFillAllCellTimeHisto)
dbba06ca 728 {
649b825d 729 Double_t time = cells->GetCellTime(absId);
dbba06ca 730 GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
649b825d 731
732 Float_t diff = (tmax-time*1.0e9);
733 fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
734 if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
735 }
736
737 }// fill cell-cluster histogram loop
738
739 }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
740
741
742 // Get vertex for photon momentum calculation and event selection
743 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 744 //GetReader()->GetVertex(v); //
649b825d 745
746 TLorentzVector mom ;
1a83b960 747 clus->GetMomentum(mom,v);
649b825d 748
749 Float_t e = mom.E();
750 Float_t pt = mom.Pt();
751 Float_t eta = mom.Eta();
752 Float_t phi = mom.Phi();
753 if(phi < 0) phi +=TMath::TwoPi();
754
755 if(GetDebug() > 0) {
756 printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
757 }
758
759 fhE ->Fill(e);
760 if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
761 if(fFillAllTH12){
762 fhPt ->Fill(pt);
763 fhPhi ->Fill(phi);
764 fhEta ->Fill(eta);
765 }
766
767 if(fFillAllTH3)
768 fhEtaPhiE->Fill(eta,phi,e);
769
770 //Cells per cluster
771 fhNCellsPerCluster ->Fill(e, nCaloCellsPerCluster);
a82b4462 772
649b825d 773 //Position
774 if(fFillAllPosHisto2){
775
776 Float_t pos[3] ;
777 clus->GetPosition(pos);
778
779 fhXE ->Fill(pos[0],e);
780 fhYE ->Fill(pos[1],e);
781 fhZE ->Fill(pos[2],e);
782 if(fFillAllTH3)
783 fhXYZ ->Fill(pos[0], pos[1],pos[2]);
784
785 fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
786 fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
787 fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
788 Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
789 fhRE ->Fill(rxyz,e);
790 fhRNCells->Fill(rxyz ,nCaloCellsPerCluster);
791 }
792
793 if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
794
795}
796
f3138ecf 797//____________________________________________________________________________
798void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters,
799 AliVCaloCells* cells)
649b825d 800{
801 // Fill clusters related histograms
802
803 TLorentzVector mom ;
804 Int_t nLabel = 0 ;
805 Int_t *labels = 0x0;
806 Int_t nCaloClusters = caloClusters->GetEntriesFast() ;
807 Int_t nCaloClustersAccepted = 0 ;
808 Int_t nCaloCellsPerCluster = 0 ;
809 Bool_t matched = kFALSE;
810 Int_t nModule =-1 ;
811
812 // Get vertex for photon momentum calculation and event selection
813 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 814 //GetReader()->GetVertex(v);
649b825d 815
816 Int_t *nClustersInModule = new Int_t[fNModules];
817 for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
818
819 if(GetDebug() > 0)
820 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
821
822 // Loop over CaloClusters
823 for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
824
1a72f6c5 825 if(GetDebug() > 0)
826 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
f3138ecf 827 iclus+1,nCaloClusters,GetReader()->GetDataType());
649b825d 828
829 AliVCluster* clus = (AliVCluster*)caloClusters->At(iclus);
830
831 // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
832 Float_t maxCellFraction = 0.;
833 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
834
835 //Cut on time of clusters
836 Double_t tof = clus->GetTOF()*1.e9;
e6fec6f5 837 if(tof < fTimeCutMin || tof > fTimeCutMax)
838 {
649b825d 839 if(GetDebug() > 0 )printf("AliAnaCalorimeterQA - Remove cluster with TOF %f\n",tof);
840 continue;
841 }
842
843 // Get cluster kinematics
1a83b960 844 clus->GetMomentum(mom,v);
649b825d 845
846 // Check only certain regions
847 Bool_t in = kTRUE;
848 if(IsFiducialCutOn()) in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
849 if(!in) continue;
850
1a83b960 851 // MC labels
649b825d 852 nLabel = clus->GetNLabels();
853 labels = clus->GetLabels();
854
1a83b960 855 // SuperModule number of cluster
856 nModule = GetModuleNumber(clus);
857
649b825d 858 // Cells per cluster
859 nCaloCellsPerCluster = clus->GetNCells();
860
861 // Cluster mathed with track?
49b5c49b 862 matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
649b825d 863
864 // Get some time averages
865 Double_t averTime[4] = {0.,0.,0.,0.};
866 CalculateAverageTime(clus, cells, averTime);
867
868 //Get time of max cell
869 Double_t tmax = cells->GetCellTime(absIdMax);
dbba06ca 870 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
649b825d 871 tmax*=1.e9;
872
1a83b960 873 // Fill histograms related to single cluster
874
875
876 // Fill some histograms before applying the exotic cell / bad map cut
877 fhNCellsPerClusterNoCut ->Fill(clus->E(), nCaloCellsPerCluster);
878 if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
879
880 fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
881
882 Float_t ampMax = cells->GetCellAmplitude(absIdMax);
dbba06ca 883 GetCaloUtils()->RecalibrateCellAmplitude(ampMax,fCalorimeter, absIdMax);
1a83b960 884
649b825d 885 //Check bad clusters if requested and rejection was not on
a82b4462 886 Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
2747966a 887
888 Float_t eCrossFrac = 0;
889 if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
649b825d 890
649b825d 891 if(!goodCluster)
1a83b960 892 {
649b825d 893 BadClusterHistograms(clus, caloClusters, cells, absIdMax,
1a83b960 894 maxCellFraction, eCrossFrac, tmax, averTime);
895 continue;
896 }
649b825d 897
898 ClusterHistograms(clus, caloClusters, cells, absIdMax,
1a83b960 899 maxCellFraction, eCrossFrac, tmax, averTime);
649b825d 900
901 nCaloClustersAccepted++;
49214ef9 902 nModule = GetModuleNumber(clus);
2747966a 903 if(nModule >=0 && nModule < fNModules)
904 {
649b825d 905 if (fCalorimeter=="EMCAL" && mom.E() > 2*fEMCALCellAmpMin) nClustersInModule[nModule]++;
906 else if(fCalorimeter=="PHOS" && mom.E() > 2*fPHOSCellAmpMin ) nClustersInModule[nModule]++;
907 }
1a83b960 908
649b825d 909 // Cluster weights
910 if(fStudyWeight) WeightHistograms(clus, cells);
911
912 // Cells in cluster position
913 if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
914
915 // Fill histograms related to single cluster, mc vs data
916 Int_t mcOK = kFALSE;
917 Int_t pdg = -1;
918 if(IsDataMC() && nLabel > 0 && labels)
919 mcOK = ClusterMCHistograms(mom, matched, labels, nLabel, pdg);
008693e5 920
649b825d 921 // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
922 if( matched && fFillAllTMHisto)
923 ClusterMatchedWithTrackHistograms(clus,mom,mcOK,pdg);
924
925 // Invariant mass
d07278cf 926 // Try to reduce background with a mild shower shape cut and no more than 1 maxima
927 // in cluster and remove low energy clusters
928 if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 &&
929 GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 &&
930 clus->GetM02() < 0.5 && clus->E() > 0.3)
a82b4462 931 InvariantMassHistograms(iclus, mom, nModule, caloClusters,cells);
649b825d 932
933 }//cluster loop
934
935 // Number of clusters histograms
936 if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
937
938 // Number of clusters per module
2747966a 939 for(Int_t imod = 0; imod < fNModules; imod++ )
940 {
649b825d 941 if(GetDebug() > 1)
942 printf("AliAnaCalorimeterQA::ClusterLoopHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]);
943 fhNClustersMod->Fill(nClustersInModule[imod],imod);
944 }
945
946 delete [] nClustersInModule;
947
948}
949
d07278cf 950//______________________________________________________________________________________________________
649b825d 951Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(const TLorentzVector mom, const Bool_t matched,
952 const Int_t * labels, const Int_t nLabels, Int_t & pdg )
953{
954
955 //Fill histograms only possible when simulation
956
2747966a 957 if(!labels || nLabels<=0)
958 {
42d47cb7 959 if(GetDebug() > 1) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Strange, labels array %p, n labels %d \n", labels,nLabels);
960 return kFALSE;
961 }
962
2747966a 963 if(GetDebug() > 1)
964 {
42d47cb7 965 printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Primaries: nlabels %d\n",nLabels);
649b825d 966 }
967
968 Float_t e = mom.E();
969 Float_t eta = mom.Eta();
970 Float_t phi = mom.Phi();
971 if(phi < 0) phi +=TMath::TwoPi();
972
973 AliAODMCParticle * aodprimary = 0x0;
974 TParticle * primary = 0x0;
975
976 //Play with the MC stack if available
977 Int_t label = labels[0];
978
2747966a 979 if(label < 0)
980 {
008693e5 981 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** bad label ***: label %d \n", label);
649b825d 982 return kFALSE;
983 }
984
985 Int_t pdg0 =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
986 Float_t vxMC= 0; Float_t vyMC = 0;
987 Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
988 Int_t charge = 0;
989
990 //Check the origin.
991 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
992
d07278cf 993 if ( GetReader()->ReadStack() &&
994 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
995 { //it MC stack and known tag
649b825d 996
d07278cf 997 if( label >= GetMCStack()->GetNtrack())
998 {
008693e5 999 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** large label ***: label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
649b825d 1000 return kFALSE;
1001 }
1002
1003 primary = GetMCStack()->Particle(label);
1004 iMother = label;
1005 pdg0 = TMath::Abs(primary->GetPdgCode());
1006 pdg = pdg0;
1007 status = primary->GetStatusCode();
1008 vxMC = primary->Vx();
1009 vyMC = primary->Vy();
1010 iParent = primary->GetFirstMother();
1011
d07278cf 1012 if(GetDebug() > 1 )
1013 {
008693e5 1014 printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Cluster most contributing mother: \n");
649b825d 1015 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1016 }
1017
1018 //Get final particle, no conversion products
d07278cf 1019 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1020 {
649b825d 1021 //Get the parent
1022 primary = GetMCStack()->Particle(iParent);
1023 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1024
1025 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Converted cluster!. Find before conversion: \n");
1026
d07278cf 1027 while((pdg == 22 || pdg == 11) && status != 1)
1028 {
008693e5 1029 Int_t iMotherOrg = iMother;
649b825d 1030 iMother = iParent;
1031 primary = GetMCStack()->Particle(iMother);
1032 status = primary->GetStatusCode();
649b825d 1033 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1034 iParent = primary->GetFirstMother();
1035
1036 // If gone too back and non stable, assign the decay photon/electron
1037 // there are other possible decays, ignore them for the moment
1038 if(pdg==111 || pdg==221)
1039 {
1040 primary = GetMCStack()->Particle(iMotherOrg);
1041 break;
1042 }
1043
1044 if( iParent < 0 )
1045 {
1046 iParent = iMother;
1047 printf("break\n");
1048 break;
1049 }
1050
649b825d 1051 if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother, primary->GetName(),status);
1052 }
008693e5 1053
d07278cf 1054 if(GetDebug() > 1 )
1055 {
649b825d 1056 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1057 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1058 }
1059
1060 }
1061
1062 //Overlapped pi0 (or eta, there will be very few), get the meson
1063 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
d07278cf 1064 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1065 {
649b825d 1066 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
008693e5 1067
d07278cf 1068 while(pdg != 111 && pdg != 221)
008693e5 1069 {
1070 //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
649b825d 1071 iMother = iParent;
1072 primary = GetMCStack()->Particle(iMother);
1073 status = primary->GetStatusCode();
649b825d 1074 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1075 iParent = primary->GetFirstMother();
1076
1077 if( iParent < 0 )break;
1078
649b825d 1079 if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg, primary->GetName(),iMother);
008693e5 1080
d07278cf 1081 if(iMother==-1)
1082 {
649b825d 1083 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1084 //break;
1085 }
1086 }
008693e5 1087
649b825d 1088 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1089 primary->GetName(),iMother);
1090 }
1091
1092 eMC = primary->Energy();
1093 ptMC = primary->Pt();
1094 phiMC = primary->Phi();
1095 etaMC = primary->Eta();
1096 pdg = TMath::Abs(primary->GetPdgCode());
1097 charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1098
1099 }
d07278cf 1100 else if( GetReader()->ReadAODMCParticles() &&
1101 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1102 {//it MC AOD and known tag
649b825d 1103 //Get the list of MC particles
1104 if(!GetReader()->GetAODMCParticles(0))
1105 AliFatal("MCParticles not available!");
1106
1107 aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1108 iMother = label;
1109 pdg0 = TMath::Abs(aodprimary->GetPdgCode());
1110 pdg = pdg0;
1111 status = aodprimary->IsPrimary();
1112 vxMC = aodprimary->Xv();
1113 vyMC = aodprimary->Yv();
1114 iParent = aodprimary->GetMother();
1115
d07278cf 1116 if(GetDebug() > 1 )
1117 {
649b825d 1118 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1119 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1120 iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1121 }
1122
1123 //Get final particle, no conversion products
008693e5 1124 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1125 {
649b825d 1126 if(GetDebug() > 1 )
1127 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1128 //Get the parent
1129 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1130 pdg = TMath::Abs(aodprimary->GetPdgCode());
d07278cf 1131 while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1132 {
008693e5 1133 Int_t iMotherOrg = iMother;
649b825d 1134 iMother = iParent;
1135 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1136 status = aodprimary->IsPrimary();
1137 iParent = aodprimary->GetMother();
1138 pdg = TMath::Abs(aodprimary->GetPdgCode());
008693e5 1139
1140 // If gone too back and non stable, assign the decay photon/electron
1141 // there are other possible decays, ignore them for the moment
1142 if(pdg==111 || pdg==221)
1143 {
1144 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMotherOrg);
1145 break;
1146 }
1147
1148 if(iParent < 0 )
1149 {
1150 iParent = iMother;
1151 break;
1152 }
1153
649b825d 1154 if(GetDebug() > 1 )
1155 printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1156 pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1157 }
1158
d07278cf 1159 if(GetDebug() > 1 )
1160 {
649b825d 1161 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1162 printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1163 iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1164 }
1165
1166 }
1167
1168 //Overlapped pi0 (or eta, there will be very few), get the meson
1169 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
2747966a 1170 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1171 {
649b825d 1172 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
008693e5 1173 while(pdg != 111 && pdg != 221)
1174 {
649b825d 1175 iMother = iParent;
1176 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1177 status = aodprimary->IsPrimary();
1178 iParent = aodprimary->GetMother();
1179 pdg = TMath::Abs(aodprimary->GetPdgCode());
008693e5 1180
1181 if(iParent < 0 ) break;
649b825d 1182
1183 if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1184
2747966a 1185 if(iMother==-1)
1186 {
649b825d 1187 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1188 //break;
1189 }
1190 }
1191
1192 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1193 aodprimary->GetName(),iMother);
1194 }
1195
1196 status = aodprimary->IsPrimary();
1197 eMC = aodprimary->E();
1198 ptMC = aodprimary->Pt();
1199 phiMC = aodprimary->Phi();
1200 etaMC = aodprimary->Eta();
1201 pdg = TMath::Abs(aodprimary->GetPdgCode());
1202 charge = aodprimary->Charge();
1203
1204 }
1205
1206 //Float_t vz = primary->Vz();
1207 Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
2747966a 1208 if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1)
1209 {
649b825d 1210 fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1211 fhEMR ->Fill(e,rVMC);
1212 }
1213
1214 //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1215 //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1216 //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1217
1218 //Overlapped pi0 (or eta, there will be very few)
2747966a 1219 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0))
1220 {
c5693f62 1221 fhRecoMCE [kmcPi0][matched] ->Fill(e,eMC);
1222 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPi0][(matched)]->Fill(eta,etaMC);
1223 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPi0][(matched)]->Fill(phi,phiMC);
1224 if(eMC > 0) fhRecoMCRatioE [kmcPi0][(matched)]->Fill(e,e/eMC);
1225 fhRecoMCDeltaE [kmcPi0][(matched)]->Fill(e,eMC-e);
1226 fhRecoMCDeltaPhi[kmcPi0][(matched)]->Fill(e,phiMC-phi);
1227 fhRecoMCDeltaEta[kmcPi0][(matched)]->Fill(e,etaMC-eta);
649b825d 1228 }//Overlapped pizero decay
2747966a 1229 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1230 {
c5693f62 1231 fhRecoMCE [kmcEta][(matched)] ->Fill(e,eMC);
1232 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcEta][(matched)]->Fill(eta,etaMC);
1233 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcEta][(matched)]->Fill(phi,phiMC);
1234 if(eMC > 0) fhRecoMCRatioE [kmcEta][(matched)]->Fill(e,e/eMC);
1235 fhRecoMCDeltaE [kmcEta][(matched)]->Fill(e,eMC-e);
1236 fhRecoMCDeltaPhi[kmcEta][(matched)]->Fill(e,phiMC-phi);
1237 fhRecoMCDeltaEta[kmcEta][(matched)]->Fill(e,etaMC-eta);
649b825d 1238 }//Overlapped eta decay
008693e5 1239 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1240 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
2747966a 1241 {
c5693f62 1242 fhRecoMCE [kmcPhoton][(matched)] ->Fill(e,eMC);
1243 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPhoton][(matched)]->Fill(eta,etaMC);
1244 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPhoton][(matched)]->Fill(phi,phiMC);
1245 if(eMC > 0) fhRecoMCRatioE [kmcPhoton][(matched)]->Fill(e,e/eMC);
008693e5 1246
c5693f62 1247 fhRecoMCDeltaE [kmcPhoton][(matched)]->Fill(e,eMC-e);
1248 fhRecoMCDeltaPhi[kmcPhoton][(matched)]->Fill(e,phiMC-phi);
1249 fhRecoMCDeltaEta[kmcPhoton][(matched)]->Fill(e,etaMC-eta);
649b825d 1250 }//photon
008693e5 1251 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1252 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
2747966a 1253 {
c5693f62 1254 fhRecoMCE [kmcElectron][(matched)] ->Fill(e,eMC);
1255 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcElectron][(matched)]->Fill(eta,etaMC);
1256 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcElectron][(matched)]->Fill(phi,phiMC);
1257 if(eMC > 0) fhRecoMCRatioE [kmcElectron][(matched)]->Fill(e,e/eMC);
1258 fhRecoMCDeltaE [kmcElectron][(matched)]->Fill(e,eMC-e);
1259 fhRecoMCDeltaPhi[kmcElectron][(matched)]->Fill(e,phiMC-phi);
1260 fhRecoMCDeltaEta[kmcElectron][(matched)]->Fill(e,etaMC-eta);
649b825d 1261 fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1262 fhEMR ->Fill(e,rVMC);
1263 }
2747966a 1264 else if(charge == 0)
1265 {
c5693f62 1266 fhRecoMCE [kmcNeHadron][(matched)] ->Fill(e,eMC);
1267 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcNeHadron][(matched)]->Fill(eta,etaMC);
1268 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcNeHadron][(matched)]->Fill(phi,phiMC);
1269 if(eMC > 0) fhRecoMCRatioE [kmcNeHadron][(matched)]->Fill(e,e/eMC);
1270 fhRecoMCDeltaE [kmcNeHadron][(matched)]->Fill(e,eMC-e);
1271 fhRecoMCDeltaPhi[kmcNeHadron][(matched)]->Fill(e,phiMC-phi);
1272 fhRecoMCDeltaEta[kmcNeHadron][(matched)]->Fill(e,etaMC-eta);
649b825d 1273 fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1274 fhHaR ->Fill(e,rVMC);
1275 }
2747966a 1276 else if(charge!=0)
1277 {
c5693f62 1278 fhRecoMCE [kmcChHadron][(matched)] ->Fill(e,eMC);
1279 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcChHadron][(matched)]->Fill(eta,etaMC);
1280 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcChHadron][(matched)]->Fill(phi,phiMC);
1281 if(eMC > 0) fhRecoMCRatioE [kmcChHadron][(matched)]->Fill(e,e/eMC);
1282 fhRecoMCDeltaE [kmcChHadron][(matched)]->Fill(e,eMC-e);
1283 fhRecoMCDeltaPhi[kmcChHadron][(matched)]->Fill(e,phiMC-phi);
1284 fhRecoMCDeltaEta[kmcChHadron][(matched)]->Fill(e,etaMC-eta);
649b825d 1285 fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1286 fhHaR ->Fill(e,rVMC);
1287 }
1288
1289 if(primary || aodprimary) return kTRUE ;
1290 else return kFALSE;
1291
1292}
1293
1294//________________________________________________________________________________________________
1295void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, TLorentzVector mom,
1296 const Bool_t okPrimary, const Int_t pdg)
1297{
1298 //Histograms for clusters matched with tracks
42d47cb7 1299
649b825d 1300 Float_t e = mom.E();
1301 Float_t pt = mom.Pt();
1302 Float_t eta = mom.Eta();
1303 Float_t phi = mom.Phi();
1304 if(phi < 0) phi +=TMath::TwoPi();
1305
2747966a 1306 if(fFillAllTH12)
1307 {
649b825d 1308 fhECharged ->Fill(e);
1309 fhPtCharged ->Fill(pt);
1310 fhPhiCharged ->Fill(phi);
1311 fhEtaCharged ->Fill(eta);
1312 }
a82b4462 1313
42d47cb7 1314 //Study the track and matched cluster if track exists.
1315
4bfeae64 1316 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
649b825d 1317
1318 if(!track) return ;
649b825d 1319
a87e069d 1320 Double_t tpt = track->Pt();
1321 Double_t tmom = track->P();
1322 Double_t dedx = track->GetTPCsignal();
1323 Int_t nITS = track->GetNcls(0);
1324 Int_t nTPC = track->GetNcls(1);
1325
1326 // Residuals
1327 Float_t deta = clus->GetTrackDz();
1328 Float_t dphi = clus->GetTrackDx();
1329 Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1330
1331 Double_t eOverP = e/tmom;
1332
d55bb5e1 1333 fh1EOverP->Fill(tpt, eOverP);
1334 if(dR < 0.02) fh1EOverPR02->Fill(tpt,eOverP);
a87e069d 1335
1336 fh2dR->Fill(e,dR);
1337 fh2MatchdEdx->Fill(tmom,dedx);
1338
1339 if(IsDataMC() && okPrimary)
1340 {
1341 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
649b825d 1342
a87e069d 1343 if(TMath::Abs(pdg) == 11)
1344 {
d55bb5e1 1345 fhMCEle1EOverP->Fill(tpt,eOverP);
a87e069d 1346 fhMCEle1dR->Fill(dR);
1347 fhMCEle2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1348 if(dR < 0.02) fhMCEle1EOverPR02->Fill(tpt,eOverP);
a87e069d 1349 }
1350 else if(charge!=0)
1351 {
d55bb5e1 1352 fhMCChHad1EOverP->Fill(tpt,eOverP);
a87e069d 1353 fhMCChHad1dR->Fill(dR);
1354 fhMCChHad2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1355 if(dR < 0.02) fhMCChHad1EOverPR02->Fill(tpt,eOverP);
a87e069d 1356 }
1357 else if(charge == 0)
1358 {
d55bb5e1 1359 fhMCNeutral1EOverP->Fill(tpt,eOverP);
a87e069d 1360 fhMCNeutral1dR->Fill(dR);
1361 fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1362 if(dR < 0.02) fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
649b825d 1363 }
a87e069d 1364 }//DataMC
1365
1366 if(dR < 0.02 && eOverP > 0.6 && eOverP< 1.2
1367 && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1368 {
1369 fh2EledEdx->Fill(tmom,dedx);
649b825d 1370 }
a87e069d 1371
649b825d 1372}
1373
1374//___________________________________
1375void AliAnaCalorimeterQA::Correlate()
1376{
1377 // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1378
1379 //Clusters
1380 TObjArray * caloClustersEMCAL = GetEMCALClusters();
1381 TObjArray * caloClustersPHOS = GetPHOSClusters();
1382
1383 Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1384 Int_t nclPHOS = caloClustersPHOS ->GetEntriesFast();
1385
1386 Float_t sumClusterEnergyEMCAL = 0;
1387 Float_t sumClusterEnergyPHOS = 0;
1388 Int_t iclus = 0;
1389 for(iclus = 0 ; iclus < caloClustersEMCAL->GetEntriesFast() ; iclus++)
1390 sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1391 for(iclus = 0 ; iclus < caloClustersPHOS->GetEntriesFast(); iclus++)
1392 sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1393
1394
1395 //Cells
1396
1397 AliVCaloCells * cellsEMCAL = GetEMCALCells();
1398 AliVCaloCells * cellsPHOS = GetPHOSCells();
1399
1400 Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1401 Int_t ncellsPHOS = cellsPHOS ->GetNumberOfCells();
1402
1403 Float_t sumCellEnergyEMCAL = 0;
1404 Float_t sumCellEnergyPHOS = 0;
1405 Int_t icell = 0;
1406 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells() ; icell++)
1407 sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1408 for(icell = 0 ; icell < cellsPHOS->GetNumberOfCells(); icell++)
1409 sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1410
1411
1412 //Fill Histograms
1413 fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1414 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1415 fhCaloCorrNCells ->Fill(ncellsEMCAL,ncellsPHOS);
1416 fhCaloCorrECells ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1417
1418 Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1419 Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1420 Int_t trM = GetTrackMultiplicity();
1421 if(fCalorimeter=="PHOS"){
1422 fhCaloV0MCorrNClusters ->Fill(v0M,nclPHOS);
1423 fhCaloV0MCorrEClusters ->Fill(v0M,sumClusterEnergyPHOS);
1424 fhCaloV0MCorrNCells ->Fill(v0M,ncellsPHOS);
1425 fhCaloV0MCorrECells ->Fill(v0M,sumCellEnergyPHOS);
1426
1427 fhCaloV0SCorrNClusters ->Fill(v0S,nclPHOS);
1428 fhCaloV0SCorrEClusters ->Fill(v0S,sumClusterEnergyPHOS);
1429 fhCaloV0SCorrNCells ->Fill(v0S,ncellsPHOS);
1430 fhCaloV0SCorrECells ->Fill(v0S,sumCellEnergyPHOS);
1431
1432 fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1433 fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1434 fhCaloTrackMCorrNCells ->Fill(trM,ncellsPHOS);
1435 fhCaloTrackMCorrECells ->Fill(trM,sumCellEnergyPHOS);
1436 }
1437 else{
1438 fhCaloV0MCorrNClusters ->Fill(v0M,nclEMCAL);
1439 fhCaloV0MCorrEClusters ->Fill(v0M,sumClusterEnergyEMCAL);
1440 fhCaloV0MCorrNCells ->Fill(v0M,ncellsEMCAL);
1441 fhCaloV0MCorrECells ->Fill(v0M,sumCellEnergyEMCAL);
1442
1443 fhCaloV0SCorrNClusters ->Fill(v0S,nclEMCAL);
1444 fhCaloV0SCorrEClusters ->Fill(v0S,sumClusterEnergyEMCAL);
1445 fhCaloV0SCorrNCells ->Fill(v0S,ncellsEMCAL);
1446 fhCaloV0SCorrECells ->Fill(v0S,sumCellEnergyEMCAL);
1447
1448 fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1449 fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1450 fhCaloTrackMCorrNCells ->Fill(trM,ncellsEMCAL);
1451 fhCaloTrackMCorrECells ->Fill(trM,sumCellEnergyEMCAL);
1452 }
1453
1454 if(GetDebug() > 0 )
1455 {
1456 printf("AliAnaCalorimeterQA::Correlate(): \n");
1457 printf("\t EMCAL: N cells %d, N clusters %d, summed E cells %f, summed E clusters %f \n",
1458 ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1459 printf("\t PHOS : N cells %d, N clusters %d, summed E cells %f, summed E clusters %f \n",
1460 ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1461 printf("\t V0 : Signal %d, Multiplicity %d, Track Multiplicity %d \n", v0S,v0M,trM);
1462 }
1463
1464}
1465
1466//__________________________________________________
1467TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1468{
1469 //Save parameters used for analysis
1470 TString parList ; //this will be list of parameters used for this analysis.
1471 const Int_t buffersize = 255;
1472 char onePar[buffersize] ;
1473
1474 snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1475 parList+=onePar ;
1476 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1477 parList+=onePar ;
1478 snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns \n",fTimeCutMin, fTimeCutMax) ;
1479 parList+=onePar ;
1480 snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1481 parList+=onePar ;
1482 //Get parameters set in base class.
1483 //parList += GetBaseParametersList() ;
1484
1485 //Get parameters set in FiducialCut class (not available yet)
1486 //parlist += GetFidCut()->GetFidCutParametersList()
1487
1488 return new TObjString(parList) ;
1489}
1490
1491//____________________________________________________
1492TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1493{
1494 // Create histograms to be saved in output file and
1495 // store them in outputContainer
1496
1497 TList * outputContainer = new TList() ;
1498 outputContainer->SetName("QAHistos") ;
1499
1500 //Histograms
745913ae 1501 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
1502 Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins(); Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax(); Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1503 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1504 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1505 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins(); Float_t massmax = GetHistogramRanges()->GetHistoMassMax(); Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
1506 Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins(); Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax(); Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
d55bb5e1 1507 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t EOverPmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t EOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
745913ae 1508 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
1509 Int_t ndRbins = GetHistogramRanges()->GetHistodRBins(); Float_t dRmax = GetHistogramRanges()->GetHistodRMax(); Float_t dRmin = GetHistogramRanges()->GetHistodRMin();
1510 Int_t ntimebins = GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1511 Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins(); Int_t nclmax = GetHistogramRanges()->GetHistoNClustersMax(); Int_t nclmin = GetHistogramRanges()->GetHistoNClustersMin();
1512 Int_t ncebins = GetHistogramRanges()->GetHistoNCellsBins(); Int_t ncemax = GetHistogramRanges()->GetHistoNCellsMax(); Int_t ncemin = GetHistogramRanges()->GetHistoNCellsMin();
1513 Int_t nceclbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nceclmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nceclmin = GetHistogramRanges()->GetHistoNClusterCellMin();
1514 Int_t nvdistbins = GetHistogramRanges()->GetHistoVertexDistBins(); Float_t vdistmax = GetHistogramRanges()->GetHistoVertexDistMax(); Float_t vdistmin = GetHistogramRanges()->GetHistoVertexDistMin();
1515 Int_t rbins = GetHistogramRanges()->GetHistoRBins(); Float_t rmax = GetHistogramRanges()->GetHistoRMax(); Float_t rmin = GetHistogramRanges()->GetHistoRMin();
1516 Int_t xbins = GetHistogramRanges()->GetHistoXBins(); Float_t xmax = GetHistogramRanges()->GetHistoXMax(); Float_t xmin = GetHistogramRanges()->GetHistoXMin();
1517 Int_t ybins = GetHistogramRanges()->GetHistoYBins(); Float_t ymax = GetHistogramRanges()->GetHistoYMax(); Float_t ymin = GetHistogramRanges()->GetHistoYMin();
1518 Int_t zbins = GetHistogramRanges()->GetHistoZBins(); Float_t zmax = GetHistogramRanges()->GetHistoZMax(); Float_t zmin = GetHistogramRanges()->GetHistoZMin();
1519 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
1520 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
1521
1522 Int_t nv0sbins = GetHistogramRanges()->GetHistoV0SignalBins(); Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax(); Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1523 Int_t nv0mbins = GetHistogramRanges()->GetHistoV0MultiplicityBins(); Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax(); Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1524 Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
649b825d 1525
1526 //EMCAL
1527 fNMaxCols = 48;
1528 fNMaxRows = 24;
1529 fNRCU = 2 ;
1530 //PHOS
1531 if(fCalorimeter=="PHOS"){
1532 fNMaxCols = 56;
1533 fNMaxRows = 64;
1534 fNRCU = 4 ;
1535 }
1536
1537 fhE = new TH1F ("hE","E reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1538 fhE->SetXTitle("E (GeV)");
1539 outputContainer->Add(fhE);
1540
1541 if(fFillAllTH12){
1542 fhPt = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1543 fhPt->SetXTitle("p_{T} (GeV/c)");
1544 outputContainer->Add(fhPt);
1545
1546 fhPhi = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1547 fhPhi->SetXTitle("#phi (rad)");
1548 outputContainer->Add(fhPhi);
1549
1550 fhEta = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1551 fhEta->SetXTitle("#eta ");
1552 outputContainer->Add(fhEta);
1553 }
1554
1555 if(fFillAllTH3){
1556 fhEtaPhiE = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1557 netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1558 fhEtaPhiE->SetXTitle("#eta ");
1559 fhEtaPhiE->SetYTitle("#phi (rad)");
1560 fhEtaPhiE->SetZTitle("E (GeV) ");
1561 outputContainer->Add(fhEtaPhiE);
1562 }
1563
1564 fhClusterTimeEnergy = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1565 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1566 fhClusterTimeEnergy->SetXTitle("E (GeV) ");
1567 fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1568 outputContainer->Add(fhClusterTimeEnergy);
1569
1570 fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1571 nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1572 fhClusterPairDiffTimeE->SetXTitle("E_{cluster} (GeV)");
1573 fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1574 outputContainer->Add(fhClusterPairDiffTimeE);
1575
1576 fhLambda0 = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1577 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1578 fhLambda0->SetXTitle("E_{cluster}");
1579 fhLambda0->SetYTitle("#lambda^{2}_{0}");
1580 outputContainer->Add(fhLambda0);
1581
1582 fhLambda1 = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1583 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1584 fhLambda1->SetXTitle("E_{cluster}");
1585 fhLambda1->SetYTitle("#lambda^{2}_{1}");
1586 outputContainer->Add(fhLambda1);
1587
1588 fhDispersion = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1589 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1590 fhDispersion->SetXTitle("E_{cluster}");
1591 fhDispersion->SetYTitle("Dispersion");
1592 outputContainer->Add(fhDispersion);
1593
1594 fhClusterMaxCellCloseCellRatio = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1595 nptbins,ptmin,ptmax, 100,0,1.);
1596 fhClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1597 fhClusterMaxCellCloseCellRatio->SetYTitle("E_{cell i}/E_{cell max}");
1598 outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1599
1600 fhClusterMaxCellCloseCellDiff = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1601 nptbins,ptmin,ptmax, 500,0,100.);
1602 fhClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1603 fhClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max}-E_{cell i} (GeV)");
1604 outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1605
1606 fhClusterMaxCellDiff = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1607 nptbins,ptmin,ptmax, 500,0,1.);
1608 fhClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1609 fhClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1610 outputContainer->Add(fhClusterMaxCellDiff);
1611
1612 fhClusterMaxCellDiffNoCut = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1613 nptbins,ptmin,ptmax, 500,0,1.);
1614 fhClusterMaxCellDiffNoCut->SetXTitle("E_{cluster} (GeV) ");
1615 fhClusterMaxCellDiffNoCut->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1616 outputContainer->Add(fhClusterMaxCellDiffNoCut);
1617
1a72f6c5 1618 fhClusterMaxCellECross = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1619 nptbins,ptmin,ptmax, 400,-1,1.);
1620 fhClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1621 fhClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1622 outputContainer->Add(fhClusterMaxCellECross);
1623
1a83b960 1624 if(fStudyBadClusters)
1625 {
649b825d 1626
1627 fhBadClusterEnergy = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1628 fhBadClusterEnergy->SetXTitle("E_{cluster} (GeV) ");
1629 outputContainer->Add(fhBadClusterEnergy);
1630
1631 fhBadClusterMaxCellCloseCellRatio = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1632 nptbins,ptmin,ptmax, 100,0,1.);
1633 fhBadClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1634 fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1635 outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1636
1637 fhBadClusterMaxCellCloseCellDiff = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1638 nptbins,ptmin,ptmax, 500,0,100);
1639 fhBadClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1640 fhBadClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max} - E_{cell i} (GeV)");
1641 outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1642
1643 fhBadClusterMaxCellDiff = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1644 nptbins,ptmin,ptmax, 500,0,1.);
1645 fhBadClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1646 fhBadClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max}) / E_{cluster}");
1647 outputContainer->Add(fhBadClusterMaxCellDiff);
1648
1649 fhBadClusterTimeEnergy = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1650 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1651 fhBadClusterTimeEnergy->SetXTitle("E_{cluster} (GeV) ");
1652 fhBadClusterTimeEnergy->SetYTitle("TOF (ns)");
1653 outputContainer->Add(fhBadClusterTimeEnergy);
1654
1655 fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1656 fhBadClusterPairDiffTimeE->SetXTitle("E_{bad cluster} (GeV)");
1657 fhBadClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1658 outputContainer->Add(fhBadClusterPairDiffTimeE);
1659
1a72f6c5 1660 fhBadClusterMaxCellECross = new TH2F ("hBadClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, bad clusters",
1661 nptbins,ptmin,ptmax, 400,-1,1.);
1662 fhBadClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1663 fhBadClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1664 outputContainer->Add(fhBadClusterMaxCellECross);
1665
e6fec6f5 1666 if(fFillAllCellTimeHisto)
1667 {
649b825d 1668 fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1669 fhBadCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
1670 fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max - i} (ns)");
1671 outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1672
1673 fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1674 fhBadClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
1675 fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
1676 outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
a82b4462 1677
649b825d 1678 fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1679 fhBadClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
1680 fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
1681 outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1682
649b825d 1683 }
1684
1685 }
1686
1687 // Cluster size in terms of cells
1688 if(fStudyClustersAsymmetry){
1689 fhDeltaIEtaDeltaIPhiE0[0] = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1690 50,0,50,50,0,50);
1691 fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
1692 fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
1693 outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
1694
1695 fhDeltaIEtaDeltaIPhiE2[0] = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1696 50,0,50,50,0,50);
1697 fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
1698 fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
1699 outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
1700
1701 fhDeltaIEtaDeltaIPhiE6[0] = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1702 50,0,50,50,0,50);
1703 fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
1704 fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
1705 outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
1706
1707 fhDeltaIA[0] = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
1708 nptbins,ptmin,ptmax,21,-1.05,1.05);
1709 fhDeltaIA[0]->SetXTitle("E_{cluster}");
1710 fhDeltaIA[0]->SetYTitle("A_{cell in cluster}");
1711 outputContainer->Add(fhDeltaIA[0]);
1712
1713 fhDeltaIAL0[0] = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
1714 ssbins,ssmin,ssmax,21,-1.05,1.05);
1715 fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
1716 fhDeltaIAL0[0]->SetYTitle("A_{cell in cluster}");
1717 outputContainer->Add(fhDeltaIAL0[0]);
1718
1719 fhDeltaIAL1[0] = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
1720 ssbins,ssmin,ssmax,21,-1.05,1.05);
1721 fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
1722 fhDeltaIAL1[0]->SetYTitle("A_{cell in cluster}");
1723 outputContainer->Add(fhDeltaIAL1[0]);
1724
1725 fhDeltaIANCells[0] = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
1726 nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1727 fhDeltaIANCells[0]->SetXTitle("N_{cell in cluster}");
1728 fhDeltaIANCells[0]->SetYTitle("A_{cell in cluster}");
1729 outputContainer->Add(fhDeltaIANCells[0]);
1730
1731
1732 fhDeltaIEtaDeltaIPhiE0[1] = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3, matched with track",
1733 50,0,50,50,0,50);
1734 fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
1735 fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
1736 outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
1737
1738 fhDeltaIEtaDeltaIPhiE2[1] = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3, matched with track",
1739 50,0,50,50,0,50);
1740 fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
1741 fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
1742 outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
1743
1744 fhDeltaIEtaDeltaIPhiE6[1] = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3, matched with track",
1745 50,0,50,50,0,50);
1746 fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
1747 fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
1748 outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
1749
1750 fhDeltaIA[1] = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
1751 nptbins,ptmin,ptmax,21,-1.05,1.05);
1752 fhDeltaIA[1]->SetXTitle("E_{cluster}");
1753 fhDeltaIA[1]->SetYTitle("A_{cell in cluster}");
1754 outputContainer->Add(fhDeltaIA[1]);
1755
1756 fhDeltaIAL0[1] = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
1757 ssbins,ssmin,ssmax,21,-1.05,1.05);
1758 fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
1759 fhDeltaIAL0[1]->SetYTitle("A_{cell in cluster}");
1760 outputContainer->Add(fhDeltaIAL0[1]);
1761
1762 fhDeltaIAL1[1] = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
1763 ssbins,ssmin,ssmax,21,-1.05,1.05);
1764 fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
1765 fhDeltaIAL1[1]->SetYTitle("A_{cell in cluster}");
1766 outputContainer->Add(fhDeltaIAL1[1]);
1767
1768 fhDeltaIANCells[1] = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
1769 nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1770 fhDeltaIANCells[1]->SetXTitle("N_{cell in cluster}");
1771 fhDeltaIANCells[1]->SetYTitle("A_{cell in cluster}");
1772 outputContainer->Add(fhDeltaIANCells[1]);
1773
1774 if(IsDataMC()){
1775 TString particle[]={"Photon","Electron","Conversion","Hadron"};
1776 for (Int_t iPart = 0; iPart < 4; iPart++) {
1777
1778 fhDeltaIAMC[iPart] = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
1779 nptbins,ptmin,ptmax,21,-1.05,1.05);
1780 fhDeltaIAMC[iPart]->SetXTitle("E_{cluster}");
1781 fhDeltaIAMC[iPart]->SetYTitle("A_{cell in cluster}");
1782 outputContainer->Add(fhDeltaIAMC[iPart]);
1783 }
1784 }
1a83b960 1785 if(fStudyBadClusters)
1786 {
1787 fhBadClusterDeltaIEtaDeltaIPhiE0 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1788 50,0,50,50,0,50);
1789 fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
1790 fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
1791 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
1792
1793 fhBadClusterDeltaIEtaDeltaIPhiE2 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1794 50,0,50,50,0,50);
1795 fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
1796 fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
1797 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
1798
1799 fhBadClusterDeltaIEtaDeltaIPhiE6 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1800 50,0,50,50,0,50);
1801 fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
1802 fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
1803 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
1804
1805 fhBadClusterDeltaIA = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
1806 nptbins,ptmin,ptmax,21,-1.05,1.05);
1807 fhBadClusterDeltaIA->SetXTitle("E_{cluster}");
1808 fhBadClusterDeltaIA->SetYTitle("A_{cell in cluster}");
1809 outputContainer->Add(fhBadClusterDeltaIA);
1810 }
649b825d 1811 }
1812
1813 if(fStudyWeight){
1814
1815 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
1816 nptbins,ptmin,ptmax, 100,0,1.);
1817 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1818 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1819 outputContainer->Add(fhECellClusterRatio);
1820
1821 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
1a72f6c5 1822 nptbins,ptmin,ptmax, 100,-10,0);
649b825d 1823 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 1824 fhECellClusterLogRatio->SetYTitle("Log(E_{cell i}/E_{cluster})");
649b825d 1825 outputContainer->Add(fhECellClusterLogRatio);
1826
1827 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
1828 nptbins,ptmin,ptmax, 100,0,1.);
1829 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1830 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1831 outputContainer->Add(fhEMaxCellClusterRatio);
1832
1833 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
1a72f6c5 1834 nptbins,ptmin,ptmax, 100,-10,0);
649b825d 1835 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 1836 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
649b825d 1837 outputContainer->Add(fhEMaxCellClusterLogRatio);
1838
1a72f6c5 1839 for(Int_t iw = 0; iw < 14; iw++){
1840 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",1+0.5*iw),
649b825d 1841 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1842 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1843 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1844 outputContainer->Add(fhLambda0ForW0[iw]);
1845
1a72f6c5 1846// fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",1+0.5*iw),
1847// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1848// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1849// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1850// outputContainer->Add(fhLambda1ForW0[iw]);
649b825d 1851
1852 if(IsDataMC()){
1853 TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
1854 for(Int_t imc = 0; imc < 5; imc++){
1855 fhLambda0ForW0MC[iw][imc] = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
1a72f6c5 1856 Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
649b825d 1857 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1858 fhLambda0ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
1859 fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
1860 outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
1861
1a72f6c5 1862// fhLambda1ForW0MC[iw][imc] = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
1863// Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
1864// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1865// fhLambda1ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
1866// fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
1867// outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
649b825d 1868 }
1869 }
1870
1871 }
1872
1873 }
1874
1875 //Track Matching
1876 if(fFillAllTMHisto){
1877 if(fFillAllTH12){
1878 fhECharged = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
1879 fhECharged->SetXTitle("E (GeV)");
1880 outputContainer->Add(fhECharged);
1881
1882 fhPtCharged = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
1883 fhPtCharged->SetXTitle("p_{T} (GeV/c)");
1884 outputContainer->Add(fhPtCharged);
1885
1886 fhPhiCharged = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
1887 fhPhiCharged->SetXTitle("#phi (rad)");
1888 outputContainer->Add(fhPhiCharged);
55c05f8c 1889
1890 fhEtaCharged = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
1891 fhEtaCharged->SetXTitle("#eta ");
1892 outputContainer->Add(fhEtaCharged);
1893 }
3748ffb5 1894 if(fFillAllTH3){
1895 fhEtaPhiECharged = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
521636d2 1896 netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
3748ffb5 1897 fhEtaPhiECharged->SetXTitle("#eta ");
1898 fhEtaPhiECharged->SetYTitle("#phi ");
1899 fhEtaPhiECharged->SetZTitle("E (GeV) ");
1900 outputContainer->Add(fhEtaPhiECharged);
1901 }
55c05f8c 1902
d55bb5e1 1903 fh1EOverP = new TH2F("h1EOverP","TRACK matches E/p",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
1904 fh1EOverP->SetYTitle("E/p");
1905 fh1EOverP->SetXTitle("p_{T} (GeV/c)");
1906 outputContainer->Add(fh1EOverP);
3bfc4732 1907
a87e069d 1908 fh2dR = new TH2F("h2dR","TRACK matches dR",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
1909 fh2dR->SetXTitle("#Delta R (rad)");
1910 fh2dR->SetXTitle("E cluster (GeV)");
1911 outputContainer->Add(fh2dR) ;
3bfc4732 1912
1913 fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
1914 fh2MatchdEdx->SetXTitle("p (GeV/c)");
1915 fh2MatchdEdx->SetYTitle("<dE/dx>");
1916 outputContainer->Add(fh2MatchdEdx);
1917
1918 fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
1919 fh2EledEdx->SetXTitle("p (GeV/c)");
1920 fh2EledEdx->SetYTitle("<dE/dx>");
1921 outputContainer->Add(fh2EledEdx) ;
1922
d55bb5e1 1923 fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches E/p, all",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
1924 fh1EOverPR02->SetYTitle("E/p");
1925 fh1EOverPR02->SetXTitle("p_{T} (GeV/c)");
1926 outputContainer->Add(fh1EOverPR02);
55c05f8c 1927 }
1928
1929 if(fFillAllPi0Histo){
9e9f04cb 1930 fhIM = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
55c05f8c 1931 fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
1932 fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
1933 outputContainer->Add(fhIM);
a6f26052 1934
55c05f8c 1935 fhAsym = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
1936 fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
1937 fhAsym->SetYTitle("Asymmetry");
1938 outputContainer->Add(fhAsym);
a6f26052 1939
a6f26052 1940 }
649b825d 1941
1a72f6c5 1942 fhNCellsPerClusterNoCut = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
35c71d5c 1943 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
715fd81f 1944 fhNCellsPerClusterNoCut->SetXTitle("E (GeV)");
1945 fhNCellsPerClusterNoCut->SetYTitle("n cells");
1946 outputContainer->Add(fhNCellsPerClusterNoCut);
649b825d 1947
1a72f6c5 1948 fhNCellsPerCluster = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
3f5990d6 1949 fhNCellsPerCluster->SetXTitle("E (GeV)");
1950 fhNCellsPerCluster->SetYTitle("n cells");
3f5990d6 1951 outputContainer->Add(fhNCellsPerCluster);
a82b4462 1952
35c71d5c 1953 fhNClusters = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
2302a644 1954 fhNClusters->SetXTitle("number of clusters");
1955 outputContainer->Add(fhNClusters);
c8fe2783 1956
55c05f8c 1957 if(fFillAllPosHisto2){
1958
1959 if(fFillAllTH3){
1960 fhXYZ = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
1961 fhXYZ->SetXTitle("x (cm)");
1962 fhXYZ->SetYTitle("y (cm)");
1963 fhXYZ->SetZTitle("z (cm) ");
1964 outputContainer->Add(fhXYZ);
1965 }
1966
35c71d5c 1967 fhXNCells = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
55c05f8c 1968 fhXNCells->SetXTitle("x (cm)");
1969 fhXNCells->SetYTitle("N cells per cluster");
1970 outputContainer->Add(fhXNCells);
1971
35c71d5c 1972 fhZNCells = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
55c05f8c 1973 fhZNCells->SetXTitle("z (cm)");
1974 fhZNCells->SetYTitle("N cells per cluster");
1975 outputContainer->Add(fhZNCells);
1976
1977 fhXE = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
1978 fhXE->SetXTitle("x (cm)");
1979 fhXE->SetYTitle("E (GeV)");
1980 outputContainer->Add(fhXE);
1981
1982 fhZE = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
1983 fhZE->SetXTitle("z (cm)");
1984 fhZE->SetYTitle("E (GeV)");
1985 outputContainer->Add(fhZE);
1986
35c71d5c 1987 fhRNCells = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
55c05f8c 1988 fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
1989 fhRNCells->SetYTitle("N cells per cluster");
1990 outputContainer->Add(fhRNCells);
1991
1992
35c71d5c 1993 fhYNCells = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
55c05f8c 1994 fhYNCells->SetXTitle("y (cm)");
1995 fhYNCells->SetYTitle("N cells per cluster");
1996 outputContainer->Add(fhYNCells);
1997
1998 fhRE = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
1999 fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2000 fhRE->SetYTitle("E (GeV)");
2001 outputContainer->Add(fhRE);
2002
2003 fhYE = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2004 fhYE->SetXTitle("y (cm)");
2005 fhYE->SetYTitle("E (GeV)");
2006 outputContainer->Add(fhYE);
2007 }
b8187de4 2008 if(fFillAllPosHisto){
521636d2 2009
a6f26052 2010 fhRCellE = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2011 fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2012 fhRCellE->SetYTitle("E (GeV)");
2013 outputContainer->Add(fhRCellE);
2014
2015 fhXCellE = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2016 fhXCellE->SetXTitle("x (cm)");
2017 fhXCellE->SetYTitle("E (GeV)");
2018 outputContainer->Add(fhXCellE);
2019
2020 fhYCellE = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2021 fhYCellE->SetXTitle("y (cm)");
2022 fhYCellE->SetYTitle("E (GeV)");
2023 outputContainer->Add(fhYCellE);
2024
2025 fhZCellE = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2026 fhZCellE->SetXTitle("z (cm)");
2027 fhZCellE->SetYTitle("E (GeV)");
2028 outputContainer->Add(fhZCellE);
2029
2030 fhXYZCell = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2031 fhXYZCell->SetXTitle("x (cm)");
2032 fhXYZCell->SetYTitle("y (cm)");
2033 fhXYZCell->SetZTitle("z (cm)");
2034 outputContainer->Add(fhXYZCell);
2035
2036
2037 Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2038 Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2039 Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2040 Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2041
35c71d5c 2042 fhDeltaCellClusterRNCells = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
a6f26052 2043 fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2044 fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
2045 outputContainer->Add(fhDeltaCellClusterRNCells);
2046
35c71d5c 2047 fhDeltaCellClusterXNCells = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
a6f26052 2048 fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
2049 fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
2050 outputContainer->Add(fhDeltaCellClusterXNCells);
2051
35c71d5c 2052 fhDeltaCellClusterYNCells = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
a6f26052 2053 fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
2054 fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2055 outputContainer->Add(fhDeltaCellClusterYNCells);
2056
35c71d5c 2057 fhDeltaCellClusterZNCells = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
a6f26052 2058 fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
2059 fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
2060 outputContainer->Add(fhDeltaCellClusterZNCells);
2061
2062 fhDeltaCellClusterRE = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2063 fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2064 fhDeltaCellClusterRE->SetYTitle("E (GeV)");
2065 outputContainer->Add(fhDeltaCellClusterRE);
2066
2067 fhDeltaCellClusterXE = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2068 fhDeltaCellClusterXE->SetXTitle("x (cm)");
2069 fhDeltaCellClusterXE->SetYTitle("E (GeV)");
2070 outputContainer->Add(fhDeltaCellClusterXE);
2071
2072 fhDeltaCellClusterYE = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2073 fhDeltaCellClusterYE->SetXTitle("y (cm)");
2074 fhDeltaCellClusterYE->SetYTitle("E (GeV)");
2075 outputContainer->Add(fhDeltaCellClusterYE);
2076
2077 fhDeltaCellClusterZE = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2078 fhDeltaCellClusterZE->SetXTitle("z (cm)");
2079 fhDeltaCellClusterZE->SetYTitle("E (GeV)");
2080 outputContainer->Add(fhDeltaCellClusterZE);
2081
2082 fhEtaPhiAmp = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2083 fhEtaPhiAmp->SetXTitle("#eta ");
2084 fhEtaPhiAmp->SetYTitle("#phi (rad)");
2085 fhEtaPhiAmp->SetZTitle("E (GeV) ");
2086 outputContainer->Add(fhEtaPhiAmp);
2087
2302a644 2088 }
c8fe2783 2089
a6f26052 2090 //Calo cells
e3300002 2091 fhNCells = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2302a644 2092 fhNCells->SetXTitle("n cells");
2093 outputContainer->Add(fhNCells);
2094
2095 fhAmplitude = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
2096 fhAmplitude->SetXTitle("Cell Energy (GeV)");
2097 outputContainer->Add(fhAmplitude);
2098
35c71d5c 2099 fhAmpId = new TH2F ("hAmpId","Cell Energy", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2302a644 2100 fhAmpId->SetXTitle("Cell Energy (GeV)");
2101 outputContainer->Add(fhAmpId);
2102
e6fec6f5 2103 if(fFillAllCellTimeHisto)
2104 {
649b825d 2105 fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
924e319f 2106 fhCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
2107 fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max-i} (ns)");
2302a644 2108 outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2109
649b825d 2110 fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
9e9f04cb 2111 fhClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
2112 fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
2113 outputContainer->Add(fhClusterMaxCellDiffAverageTime);
a82b4462 2114
649b825d 2115 fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2116 fhClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
2117 fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
2118 outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
35c71d5c 2119
924e319f 2120 fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
35c71d5c 2121 fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2302a644 2122 fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2123 outputContainer->Add(fhCellIdCellLargeTimeSpread);
521636d2 2124
2302a644 2125 fhTime = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
2126 fhTime->SetXTitle("Cell Time (ns)");
2127 outputContainer->Add(fhTime);
2128
1a72f6c5 2129 fhTimeVz = new TH2F ("hTimeVz","Cell Time vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2130 fhTimeVz->SetXTitle("|v_{z}| (cm)");
2131 fhTimeVz->SetYTitle("Cell Time (ns)");
2132 outputContainer->Add(fhTimeVz);
2133
35c71d5c 2134 fhTimeId = new TH2F ("hTimeId","Cell Time vs Absolute Id",
2135 ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2302a644 2136 fhTimeId->SetXTitle("Cell Time (ns)");
2137 fhTimeId->SetYTitle("Cell Absolute Id");
2138 outputContainer->Add(fhTimeId);
2139
2140 fhTimeAmp = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2141 fhTimeAmp->SetYTitle("Cell Time (ns)");
2142 fhTimeAmp->SetXTitle("Cell Energy (GeV)");
2143 outputContainer->Add(fhTimeAmp);
2144
2302a644 2145 }
06f1b12a 2146
2147 fhCellECross = new TH2F ("hCellECross","1 - Energy in cross around cell / cell energy",
2148 nptbins,ptmin,ptmax, 400,-1,1.);
2149 fhCellECross->SetXTitle("E_{cell} (GeV) ");
2150 fhCellECross->SetYTitle("1- E_{cross}/E_{cell}");
2151 outputContainer->Add(fhCellECross);
2152
1a72f6c5 2153
798a9b04 2154 if(fCorrelate){
2155 //PHOS vs EMCAL
35c71d5c 2156 fhCaloCorrNClusters = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2302a644 2157 fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2158 fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2159 outputContainer->Add(fhCaloCorrNClusters);
2160
798a9b04 2161 fhCaloCorrEClusters = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2302a644 2162 fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
2163 fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
2164 outputContainer->Add(fhCaloCorrEClusters);
2165
35c71d5c 2166 fhCaloCorrNCells = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
649b825d 2167 fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2168 fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2169 outputContainer->Add(fhCaloCorrNCells);
2302a644 2170
649b825d 2171 fhCaloCorrECells = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins*2,ptmin,ptmax*2);
2172 fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
2173 fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
2174 outputContainer->Add(fhCaloCorrECells);
2302a644 2175
649b825d 2176 //Calorimeter VS V0 signal
2177 fhCaloV0SCorrNClusters = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2178 fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2179 fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2180 outputContainer->Add(fhCaloV0SCorrNClusters);
2302a644 2181
649b825d 2182 fhCaloV0SCorrEClusters = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2183 fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2184 fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2185 outputContainer->Add(fhCaloV0SCorrEClusters);
2302a644 2186
649b825d 2187 fhCaloV0SCorrNCells = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2188 fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2189 fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2190 outputContainer->Add(fhCaloV0SCorrNCells);
3bfc4732 2191
649b825d 2192 fhCaloV0SCorrECells = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2193 fhCaloV0SCorrECells->SetXTitle("V0 signal");
2194 fhCaloV0SCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2195 outputContainer->Add(fhCaloV0SCorrECells);
3bfc4732 2196
649b825d 2197 //Calorimeter VS V0 multiplicity
2198 fhCaloV0MCorrNClusters = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2199 fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2200 fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2201 outputContainer->Add(fhCaloV0MCorrNClusters);
3bfc4732 2202
649b825d 2203 fhCaloV0MCorrEClusters = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2204 fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2205 fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2206 outputContainer->Add(fhCaloV0MCorrEClusters);
3bfc4732 2207
649b825d 2208 fhCaloV0MCorrNCells = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2209 fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2210 fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2211 outputContainer->Add(fhCaloV0MCorrNCells);
3bfc4732 2212
649b825d 2213 fhCaloV0MCorrECells = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2214 fhCaloV0MCorrECells->SetXTitle("V0 signal");
2215 fhCaloV0MCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2216 outputContainer->Add(fhCaloV0MCorrECells);
3bfc4732 2217
649b825d 2218 //Calorimeter VS Track multiplicity
2219 fhCaloTrackMCorrNClusters = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2220 fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2221 fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2222 outputContainer->Add(fhCaloTrackMCorrNClusters);
3bfc4732 2223
649b825d 2224 fhCaloTrackMCorrEClusters = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2225 fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2226 fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2227 outputContainer->Add(fhCaloTrackMCorrEClusters);
3bfc4732 2228
649b825d 2229 fhCaloTrackMCorrNCells = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2230 fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2231 fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2232 outputContainer->Add(fhCaloTrackMCorrNCells);
3bfc4732 2233
649b825d 2234 fhCaloTrackMCorrECells = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2235 fhCaloTrackMCorrECells->SetXTitle("# tracks");
2236 fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2237 outputContainer->Add(fhCaloTrackMCorrECells);
3bfc4732 2238
3bfc4732 2239
649b825d 2240 }//correlate calorimeters
c8fe2783 2241
649b825d 2242 //Module histograms
9725fd2a 2243
649b825d 2244 fhEMod = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2245 fhEMod->SetXTitle("E (GeV)");
2246 fhEMod->SetYTitle("Module");
2247 outputContainer->Add(fhEMod);
c8fe2783 2248
649b825d 2249 fhAmpMod = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2250 fhAmpMod->SetXTitle("E (GeV)");
2251 fhAmpMod->SetYTitle("Module");
2252 outputContainer->Add(fhAmpMod);
3bfc4732 2253
e6fec6f5 2254 if(fFillAllCellTimeHisto)
2255 {
0fb69ade 2256 fhTimeMod = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2257 fhTimeMod->SetXTitle("t (ns)");
2258 fhTimeMod->SetYTitle("Module");
2259 outputContainer->Add(fhTimeMod);
2260 }
3bfc4732 2261
e3300002 2262 fhNClustersMod = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
649b825d 2263 fhNClustersMod->SetXTitle("number of clusters");
2264 fhNClustersMod->SetYTitle("Module");
2265 outputContainer->Add(fhNClustersMod);
2302a644 2266
e3300002 2267 fhNCellsMod = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
649b825d 2268 fhNCellsMod->SetXTitle("n cells");
2269 fhNCellsMod->SetYTitle("Module");
2270 outputContainer->Add(fhNCellsMod);
17708df9 2271
649b825d 2272 Int_t colmaxs = fNMaxCols;
2273 Int_t rowmaxs = fNMaxRows;
e6fec6f5 2274 if(fCalorimeter=="EMCAL")
2275 {
649b825d 2276 colmaxs=2*fNMaxCols;
2277 rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2278 }
e6fec6f5 2279 else
2280 {
649b825d 2281 rowmaxs=fNModules*fNMaxRows;
2302a644 2282 }
3bfc4732 2283
649b825d 2284 fhGridCells = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2285 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2286 fhGridCells->SetYTitle("row (phi direction)");
2287 fhGridCells->SetXTitle("column (eta direction)");
2288 outputContainer->Add(fhGridCells);
3bfc4732 2289
649b825d 2290 fhGridCellsE = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2291 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2292 fhGridCellsE->SetYTitle("row (phi direction)");
2293 fhGridCellsE->SetXTitle("column (eta direction)");
2294 outputContainer->Add(fhGridCellsE);
3bfc4732 2295
e6fec6f5 2296 if(fFillAllCellTimeHisto)
2297 {
0fb69ade 2298 fhGridCellsTime = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2299 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2300 fhGridCellsTime->SetYTitle("row (phi direction)");
2301 fhGridCellsTime->SetXTitle("column (eta direction)");
2302 outputContainer->Add(fhGridCellsTime);
2303 }
3bfc4732 2304
e6fec6f5 2305 fhNCellsPerClusterMod = new TH2F*[fNModules];
649b825d 2306 fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
e6fec6f5 2307 fhIMMod = new TH2F*[fNModules];
2308 if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
0fb69ade 2309
e6fec6f5 2310 for(Int_t imod = 0; imod < fNModules; imod++)
2311 {
649b825d 2312 fhNCellsPerClusterMod[imod] = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2313 Form("# cells per cluster vs cluster energy in Module %d",imod),
2314 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2315 fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
2316 fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
2317 outputContainer->Add(fhNCellsPerClusterMod[imod]);
3bfc4732 2318
649b825d 2319 fhNCellsPerClusterModNoCut[imod] = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2320 Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2321 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2322 fhNCellsPerClusterModNoCut[imod]->SetXTitle("E (GeV)");
2323 fhNCellsPerClusterModNoCut[imod]->SetYTitle("n cells");
2324 outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
3bfc4732 2325
e6fec6f5 2326 if(fFillAllCellTimeHisto)
2327 {
2328 for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2329 {
649b825d 2330 fhTimeAmpPerRCU[imod*fNRCU+ircu] = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2331 Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu),
2332 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2333 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
2334 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
2335 outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
3bfc4732 2336
649b825d 2337 }
3bfc4732 2338 }
2339
e6fec6f5 2340 if(fFillAllPi0Histo)
2341 {
649b825d 2342 fhIMMod[imod] = new TH2F (Form("hIM_Mod%d",imod),
2343 Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2344 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2345 fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
2346 fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2347 outputContainer->Add(fhIMMod[imod]);
3bfc4732 2348
649b825d 2349 }
2350 }
2351
a82b4462 2352 // Monte Carlo Histograms
649b825d 2353
2354 TString particleName[] = { "Photon", "Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2355
2356 if(IsDataMC()){
2357 for(Int_t iPart = 0; iPart < 6; iPart++){
a82b4462 2358
649b825d 2359 for(Int_t iCh = 0; iCh < 2; iCh++){
3bfc4732 2360
649b825d 2361 fhRecoMCRatioE[iPart][iCh] = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2362 Form("Reco/Gen E, %s, Matched %d",particleName[iPart].Data(),iCh),
2363 nptbins, ptmin, ptmax, 200,0,2);
2364 fhRecoMCRatioE[iPart][iCh]->SetXTitle("E_{reconstructed}/E_{generated}");
2365 outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
3bfc4732 2366
3bfc4732 2367
649b825d 2368 fhRecoMCDeltaE[iPart][iCh] = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
2369 Form("MC - Reco E, %s, Matched %d",particleName[iPart].Data(),iCh),
2370 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
2371 fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#Delta E (GeV)");
2372 outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
3bfc4732 2373
649b825d 2374 fhRecoMCDeltaPhi[iPart][iCh] = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2375 Form("MC - Reco #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2376 nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
2377 fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#Delta #phi (rad)");
2378 outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
3bfc4732 2379
649b825d 2380 fhRecoMCDeltaEta[iPart][iCh] = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2381 Form("MC- Reco #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2382 nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
2383 fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#Delta #eta ");
2384 outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
e1e62b89 2385
649b825d 2386 fhRecoMCE[iPart][iCh] = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2387 Form("E distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2388 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2389 fhRecoMCE[iPart][iCh]->SetXTitle("E_{rec} (GeV)");
2390 fhRecoMCE[iPart][iCh]->SetYTitle("E_{gen} (GeV)");
2391 outputContainer->Add(fhRecoMCE[iPart][iCh]);
3bfc4732 2392
649b825d 2393 fhRecoMCPhi[iPart][iCh] = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2394 Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2395 nphibins,phimin,phimax, nphibins,phimin,phimax);
2396 fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{rec} (rad)");
2397 fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{gen} (rad)");
2398 outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
e1e62b89 2399
649b825d 2400 fhRecoMCEta[iPart][iCh] = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2401 Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2402 netabins,etamin,etamax,netabins,etamin,etamax);
2403 fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{rec} ");
2404 fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{gen} ");
2405 outputContainer->Add(fhRecoMCEta[iPart][iCh]);
3bfc4732 2406 }
649b825d 2407 }
c8fe2783 2408
649b825d 2409 //Pure MC
2410 for(Int_t iPart = 0; iPart < 4; iPart++){
2411 fhGenMCE[iPart] = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2412 Form("p_{T} of generated %s",particleName[iPart].Data()),
2413 nptbins,ptmin,ptmax);
2414 fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2415 Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2416 netabins,etamin,etamax,nphibins,phimin,phimax);
2417
2418 fhGenMCE[iPart] ->SetXTitle("p_{T} (GeV/c)");
2419 fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2420 fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
3bfc4732 2421
649b825d 2422 outputContainer->Add(fhGenMCE[iPart]);
2423 outputContainer->Add(fhGenMCEtaPhi[iPart]);
521636d2 2424
521636d2 2425
649b825d 2426 fhGenMCAccE[iPart] = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2427 Form("p_{T} of generated %s",particleName[iPart].Data()),
2428 nptbins,ptmin,ptmax);
2429 fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2430 Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2431 netabins,etamin,etamax,nphibins,phimin,phimax);
2432
2433 fhGenMCAccE[iPart] ->SetXTitle("p_{T} (GeV/c)");
2434 fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2435 fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2436
2437 outputContainer->Add(fhGenMCAccE[iPart]);
2438 outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2439
2440 }
f5036bcb 2441
649b825d 2442 //Vertex of generated particles
2443
2444 fhEMVxyz = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2445 fhEMVxyz->SetXTitle("v_{x}");
2446 fhEMVxyz->SetYTitle("v_{y}");
2447 //fhEMVxyz->SetZTitle("v_{z}");
2448 outputContainer->Add(fhEMVxyz);
2449
2450 fhHaVxyz = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2451 fhHaVxyz->SetXTitle("v_{x}");
2452 fhHaVxyz->SetYTitle("v_{y}");
2453 //fhHaVxyz->SetZTitle("v_{z}");
2454 outputContainer->Add(fhHaVxyz);
2455
2456 fhEMR = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2457 fhEMR->SetXTitle("E (GeV)");
2458 fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2459 outputContainer->Add(fhEMR);
2460
2461 fhHaR = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2462 fhHaR->SetXTitle("E (GeV)");
2463 fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2464 outputContainer->Add(fhHaR);
2465
2466
2467 //Track Matching
2468
d55bb5e1 2469 fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2470 fhMCEle1EOverP->SetYTitle("E/p");
2471 fhMCEle1EOverP->SetXTitle("p_{T} (GeV/c)");
2472 outputContainer->Add(fhMCEle1EOverP);
649b825d 2473
2474 fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2475 fhMCEle1dR->SetXTitle("#Delta R (rad)");
2476 outputContainer->Add(fhMCEle1dR) ;
2477
2478 fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2479 fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
2480 fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
2481 outputContainer->Add(fhMCEle2MatchdEdx);
2482
d55bb5e1 2483 fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2484 fhMCChHad1EOverP->SetYTitle("E/p");
2485 fhMCChHad1EOverP->SetXTitle("p_{T} (GeV/c)");
2486 outputContainer->Add(fhMCChHad1EOverP);
649b825d 2487
2488 fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2489 fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2490 outputContainer->Add(fhMCChHad1dR) ;
2491
2492 fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2493 fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
2494 fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
2495 outputContainer->Add(fhMCChHad2MatchdEdx);
2496
d55bb5e1 2497 fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2498 fhMCNeutral1EOverP->SetYTitle("E/p");
2499 fhMCNeutral1EOverP->SetXTitle("p_{T} (GeV/c)");
2500 outputContainer->Add(fhMCNeutral1EOverP);
649b825d 2501
2502 fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2503 fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
2504 outputContainer->Add(fhMCNeutral1dR) ;
2505
2506 fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2507 fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
2508 fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
2509 outputContainer->Add(fhMCNeutral2MatchdEdx);
2510
d55bb5e1 2511 fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2512 fhMCEle1EOverPR02->SetYTitle("E/p");
2513 fhMCEle1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2514 outputContainer->Add(fhMCEle1EOverPR02);
649b825d 2515
d55bb5e1 2516 fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2517 fhMCChHad1EOverPR02->SetYTitle("E/p");
2518 fhMCChHad1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2519 outputContainer->Add(fhMCChHad1EOverPR02);
f5036bcb 2520
d55bb5e1 2521 fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2522 fhMCNeutral1EOverPR02->SetYTitle("E/p");
2523 fhMCNeutral1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2524 outputContainer->Add(fhMCNeutral1EOverPR02);
521636d2 2525 }
f5036bcb 2526
649b825d 2527 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
2528 // printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
521636d2 2529
649b825d 2530 return outputContainer;
902aa95c 2531}
2532
1a72f6c5 2533//_____________________________________________________________________________________________
2534Float_t AliAnaCalorimeterQA::GetECross(const Int_t absID, AliVCaloCells* cells)
2535{
2536 // Get energy in cross axis around maximum cell, for EMCAL only
2537
06f1b12a 2538 Int_t icol =-1, irow=-1,iRCU = -1;
2539 Int_t imod = GetModuleNumberCellIndexes(absID, fCalorimeter, icol, irow, iRCU);
e3300002 2540
57d8227a 2541 if(fCalorimeter=="EMCAL")
2542 {
06f1b12a 2543 //Get close cells index, energy and time, not in corners
e3300002 2544
2545 Int_t absID1 = -1;
2546 Int_t absID2 = -1;
2547
2548 if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
2549 if( irow > 0 ) absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
2747966a 2550
2551 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2552 Int_t absID3 = -1;
2553 Int_t absID4 = -1;
2554
2555 if ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
2556 {
e3300002 2557 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
2558 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol-1);
2747966a 2559 }
2560 else if( icol == 0 && imod%2 )
2561 {
e3300002 2562 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol+1);
2563 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
2747966a 2564 }
2565 else
2566 {
e3300002 2567 if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
2568 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
2569 if( icol > 0 )
2570 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
2747966a 2571 }
06f1b12a 2572
2573 //Recalibrate cell energy if needed
2574 //Float_t ecell = cells->GetCellAmplitude(absID);
dbba06ca 2575 //GetCaloUtils()->RecalibrateCellAmplitude(ecell,fCalorimeter, absID);
06f1b12a 2576 Double_t tcell = cells->GetCellTime(absID);
dbba06ca 2577 GetCaloUtils()->RecalibrateCellTime(tcell, fCalorimeter, absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2578
2579 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2580 Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
2581
2747966a 2582 if(absID1 >0 )
2583 {
06f1b12a 2584 ecell1 = cells->GetCellAmplitude(absID1);
dbba06ca 2585 GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absID1);
06f1b12a 2586 tcell1 = cells->GetCellTime(absID1);
dbba06ca 2587 GetCaloUtils()->RecalibrateCellTime (tcell1, fCalorimeter, absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2588 }
2747966a 2589 if(absID2 >0 )
2590 {
06f1b12a 2591 ecell2 = cells->GetCellAmplitude(absID2);
dbba06ca 2592 GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absID2);
06f1b12a 2593 tcell2 = cells->GetCellTime(absID2);
dbba06ca 2594 GetCaloUtils()->RecalibrateCellTime (tcell2, fCalorimeter, absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2595 }
2747966a 2596 if(absID3 >0 )
2597 {
06f1b12a 2598 ecell3 = cells->GetCellAmplitude(absID3);
dbba06ca 2599 GetCaloUtils()->RecalibrateCellAmplitude(ecell3, fCalorimeter, absID3);
06f1b12a 2600 tcell3 = cells->GetCellTime(absID3);
dbba06ca 2601 GetCaloUtils()->RecalibrateCellTime (tcell3, fCalorimeter, absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2602 }
2747966a 2603 if(absID4 >0 )
2604 {
06f1b12a 2605 ecell4 = cells->GetCellAmplitude(absID4);
dbba06ca 2606 GetCaloUtils()->RecalibrateCellAmplitude(ecell4, fCalorimeter, absID4);
06f1b12a 2607 tcell4 = cells->GetCellTime(absID4);
dbba06ca 2608 GetCaloUtils()->RecalibrateCellTime (tcell4, fCalorimeter, absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2609 }
2610
2611 if(TMath::Abs(tcell-tcell1)*1.e9 > 50) ecell1 = 0 ;
2612 if(TMath::Abs(tcell-tcell2)*1.e9 > 50) ecell2 = 0 ;
2613 if(TMath::Abs(tcell-tcell3)*1.e9 > 50) ecell3 = 0 ;
2614 if(TMath::Abs(tcell-tcell4)*1.e9 > 50) ecell4 = 0 ;
2615
2616 return ecell1+ecell2+ecell3+ecell4;
1a72f6c5 2617 }
57d8227a 2618 else //PHOS
2619 {
06f1b12a 2620
2621 Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
2622
2623 Int_t relId1[] = { imod+1, 0, irow+1, icol };
2624 Int_t relId2[] = { imod+1, 0, irow-1, icol };
2625 Int_t relId3[] = { imod+1, 0, irow , icol+1 };
2626 Int_t relId4[] = { imod+1, 0, irow , icol-1 };
2627
2628 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
2629 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
2630 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
2631 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
2632
2633 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2634
2635 if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
2636 if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
2637 if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
2638 if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
2639
2640 return ecell1+ecell2+ecell3+ecell4;
2641
1a72f6c5 2642 }
2643
1a72f6c5 2644}
2645
c5693f62 2646//__________________________________________________________________________________________________
649b825d 2647void AliAnaCalorimeterQA::InvariantMassHistograms(const Int_t iclus, const TLorentzVector mom,
c5693f62 2648 const Int_t nModule, const TObjArray* caloClusters,
a82b4462 2649 AliVCaloCells * cells)
649b825d 2650{
2651 // Fill Invariant mass histograms
c8fe2783 2652
649b825d 2653 if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
3748ffb5 2654
649b825d 2655 //Get vertex for photon momentum calculation and event selection
2656 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 2657 //GetReader()->GetVertex(v);
a6f26052 2658
649b825d 2659 Int_t nModule2 = -1;
2660 TLorentzVector mom2 ;
2661 Int_t nCaloClusters = caloClusters->GetEntriesFast();
a6f26052 2662
d07278cf 2663 for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
2664 {
649b825d 2665 AliVCluster* clus2 = (AliVCluster*)caloClusters->At(jclus);
a82b4462 2666
2667 Float_t maxCellFraction = 0.;
2668 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction);
55c05f8c 2669
d07278cf 2670 // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
2671 // in cluster and remove low energy clusters
2672 if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
2673 GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
2674 clus2->GetM02() > 0.5 || clus2->E() < 0.3) continue;
c8fe2783 2675
649b825d 2676 //Get cluster kinematics
2677 clus2->GetMomentum(mom2,v);
c8fe2783 2678
649b825d 2679 //Check only certain regions
2680 Bool_t in2 = kTRUE;
2681 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
2682 if(!in2) continue;
2302a644 2683
649b825d 2684 //Get module of cluster
2685 nModule2 = GetModuleNumber(clus2);
c8fe2783 2686
649b825d 2687 //Fill histograms
c8fe2783 2688
649b825d 2689 //All modules
2690 fhIM ->Fill((mom+mom2).Pt(),(mom+mom2).M());
49214ef9 2691
649b825d 2692 //Single module
d07278cf 2693 if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
649b825d 2694 fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
d07278cf 2695
c8fe2783 2696
649b825d 2697 //Asymetry histograms
2698 fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
2699
2700 }// 2nd cluster loop
2701
2702}
2703
2704//______________________________
2705void AliAnaCalorimeterQA::Init()
2706{
2707 //Check if the data or settings are ok
c8fe2783 2708
649b825d 2709 if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL")
2710 AliFatal(Form("Wrong calorimeter name <%s>", fCalorimeter.Data()));
521636d2 2711
649b825d 2712 if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
2713 AliFatal("Analysis of reconstructed data, MC reader not aplicable");
2714
2715}
3bfc4732 2716
649b825d 2717//________________________________________
2718void AliAnaCalorimeterQA::InitParameters()
2719{
2720 //Initialize the parameters of the analysis.
2721 AddToHistogramsName("AnaCaloQA_");
2722
2723 fCalorimeter = "EMCAL"; //or PHOS
2724 fNModules = 12; // set maximum to maximum number of EMCAL modules
2725 fNRCU = 2; // set maximum number of RCU in EMCAL per SM
e6fec6f5 2726 fTimeCutMin = -9999999;
2727 fTimeCutMax = 9999999;
649b825d 2728 fEMCALCellAmpMin = 0.2;
2729 fPHOSCellAmpMin = 0.2;
c8fe2783 2730
649b825d 2731}
c8fe2783 2732
a82b4462 2733//___________________________________________________________________________________
2734Bool_t AliAnaCalorimeterQA::IsGoodCluster(const Int_t absIdMax, AliVCaloCells* cells)
649b825d 2735{
2736 //Identify cluster as exotic or not
2737
06f1b12a 2738 if(!fStudyBadClusters) return kTRUE;
a82b4462 2739
2747966a 2740 if(fCalorimeter=="EMCAL")
2741 {
06f1b12a 2742 if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
2747966a 2743 {
06f1b12a 2744 return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
2747966a 2745 }
2746 else
2747 {
06f1b12a 2748 return kTRUE;
2747966a 2749 }
649b825d 2750 }
06f1b12a 2751 else // PHOS
2752 {
1a83b960 2753 Float_t ampMax = cells->GetCellAmplitude(absIdMax);
dbba06ca 2754 GetCaloUtils()->RecalibrateCellAmplitude(ampMax, fCalorimeter, absIdMax);
2747966a 2755
2756 if(ampMax < 0.01) return kFALSE;
2757
1a83b960 2758 if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
2759 else return kTRUE;
06f1b12a 2760 }
2761
649b825d 2762}
17708df9 2763
a82b4462 2764//_________________________________________________________
649b825d 2765void AliAnaCalorimeterQA::Print(const Option_t * opt) const
2766{
2767 //Print some relevant parameters set for the analysis
2768 if(! opt)
2769 return;
521636d2 2770
649b825d 2771 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2772 AliAnaCaloTrackCorrBaseClass::Print(" ");
2302a644 2773
649b825d 2774 printf("Select Calorimeter %s \n",fCalorimeter.Data());
2775 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2776 printf("EMCAL Min Amplitude : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
2777 printf("PHOS Min Amplitude : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
c8fe2783 2778
649b825d 2779}
2780
649b825d 2781//_____________________________________________________
2782void AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
2783{
2784 //Fill Calorimeter QA histograms
c8fe2783 2785
649b825d 2786 //Play with the MC stack if available
2787 if(IsDataMC()) MCHistograms();
798a9b04 2788
649b825d 2789 //Get List with CaloClusters
2790 TObjArray * caloClusters = NULL;
2791 if (fCalorimeter == "PHOS") caloClusters = GetPHOSClusters();
2792 else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
2793 else
2794 AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
798a9b04 2795
649b825d 2796 // Do not do anything if there are no clusters
2797 if(caloClusters->GetEntriesFast() == 0) return;
521636d2 2798
649b825d 2799 //Get List with CaloCells
2800 AliVCaloCells * cells = 0x0;
2801 if(fCalorimeter == "PHOS") cells = GetPHOSCells();
2802 else cells = GetEMCALCells();
798a9b04 2803
649b825d 2804 if(!caloClusters || !cells) {
2805 AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
2806 return; // trick coverity
c8fe2783 2807 }
649b825d 2808
1a72f6c5 2809 //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
2810
649b825d 2811 // Correlate Calorimeters and V0 and track Multiplicity
2812 if(fCorrelate) Correlate();
2813
2814 // Clusters
2815 ClusterLoopHistograms(caloClusters,cells);
2816
2817 // Cells
2818 CellHistograms(cells);
2819
2820 if(GetDebug() > 0)
2821 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
2822
a0bb4dc0 2823}
2824
649b825d 2825//______________________________________
2826void AliAnaCalorimeterQA::MCHistograms()
2827{
2828 //Get the MC arrays and do some checks before filling MC histograms
9e9f04cb 2829
649b825d 2830 TLorentzVector mom ;
2831
2832 if(GetReader()->ReadStack()){
9e9f04cb 2833
649b825d 2834 if(!GetMCStack())
2835 AliFatal("Stack not available, is the MC handler called?\n");
9e9f04cb 2836
649b825d 2837 //Fill some pure MC histograms, only primaries.
2747966a 2838 for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++)
2839 {//Only primary particles, for all MC transport put GetNtrack()
649b825d 2840 TParticle *primary = GetMCStack()->Particle(i) ;
9e9f04cb 2841
649b825d 2842 if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
2843 primary->Momentum(mom);
2844 MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
2845 } //primary loop
2846 }
2847 else if(GetReader()->ReadAODMCParticles()){
2848
2849 if(!GetReader()->GetAODMCParticles(0))
2850 AliFatal("AODMCParticles not available!");
2851
2852 //Fill some pure MC histograms, only primaries.
2747966a 2853 for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++)
2854 {
649b825d 2855 AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
9e9f04cb 2856
649b825d 2857 if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
9e9f04cb 2858
649b825d 2859 mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
2860 MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
2861 } //primary loop
2862
2863 }
2864}
17708df9 2865
649b825d 2866//_______________________________________________________________________________
2867void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg)
2868{
a6f26052 2869 //Fill pure monte carlo related histograms
4865e325 2870
2302a644 2871 Float_t eMC = mom.E();
2302a644 2872 Float_t phiMC = mom.Phi();
2873 if(phiMC < 0)
2874 phiMC += TMath::TwoPi();
2875 Float_t etaMC = mom.Eta();
2876
2877 if (TMath::Abs(etaMC) > 1) return;
2878
35c71d5c 2879 Bool_t in = kFALSE;
2880
2881 //Rough stimate of acceptance for pi0, Eta and electrons
2747966a 2882 if(fCalorimeter == "PHOS")
2883 {
35c71d5c 2884 if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
2885 in = kTRUE ;
2886 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
2887
2888 }
2747966a 2889 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
2890 {
2891 if(GetEMCALGeometry())
2892 {
35c71d5c 2893 Int_t absID=0;
2894 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID);
2895
2896 if( absID >= 0)
2897 in = kTRUE;
2898
2899 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in);
2900 }
2747966a 2901 else
2902 {
35c71d5c 2903 if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
2904 in = kTRUE ;
2905 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
2906 }
2907 }
2302a644 2908
2747966a 2909 if (pdg==22)
2910 {
c5693f62 2911 fhGenMCE[kmcPhoton] ->Fill(eMC);
2912 if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
2747966a 2913 if(in)
2914 {
c5693f62 2915 fhGenMCAccE[kmcPhoton] ->Fill(eMC);
2916 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
2302a644 2917 }
2918 }
2747966a 2919 else if (pdg==111)
2920 {
c5693f62 2921 fhGenMCE[kmcPi0] ->Fill(eMC);
2922 if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
2747966a 2923 if(in)
2924 {
c5693f62 2925 fhGenMCAccE[kmcPi0] ->Fill(eMC);
2926 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
2302a644 2927 }
2928 }
2747966a 2929 else if (pdg==221)
2930 {
c5693f62 2931 fhGenMCE[kmcEta] ->Fill(eMC);
2932 if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC);
2747966a 2933 if(in)
2934 {
c5693f62 2935 fhGenMCAccE[kmcEta] ->Fill(eMC);
2936 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC);
35c71d5c 2937 }
2302a644 2938 }
2747966a 2939 else if (TMath::Abs(pdg)==11)
2940 {
c5693f62 2941 fhGenMCE[kmcElectron] ->Fill(eMC);
2942 if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
2747966a 2943 if(in)
2944 {
c5693f62 2945 fhGenMCAccE[kmcElectron] ->Fill(eMC);
2946 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
35c71d5c 2947 }
2302a644 2948 }
902aa95c 2949}
c8fe2783 2950
649b825d 2951//_________________________________________________________________________________
2952void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
2953{
2954 // Calculate weights
2955
2956 // First recalculate energy in case non linearity was applied
2957 Float_t energy = 0;
2958 Float_t ampMax = 0;
2747966a 2959 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
2960 {
649b825d 2961 Int_t id = clus->GetCellsAbsId()[ipos];
2962
2963 //Recalibrate cell energy if needed
2964 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 2965 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 2966
2967 energy += amp;
2968
2969 if(amp> ampMax)
2970 ampMax = amp;
2971
2972 } // energy loop
2973
2747966a 2974 if(energy <=0 )
2975 {
649b825d 2976 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
2977 return;
2978 }
2979
2980 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
2981 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
2982
2983 //Get the ratio and log ratio to all cells in cluster
2747966a 2984 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
2985 {
649b825d 2986 Int_t id = clus->GetCellsAbsId()[ipos];
2987
2988 //Recalibrate cell energy if needed
2989 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 2990 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 2991
2992 fhECellClusterRatio ->Fill(energy,amp/energy);
2993 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
2994 }
2995
2996 //Recalculate shower shape for different W0
2747966a 2997 if(fCalorimeter=="EMCAL")
2998 {
649b825d 2999 Float_t l0org = clus->GetM02();
3000 Float_t l1org = clus->GetM20();
3001 Float_t dorg = clus->GetDispersion();
3002
1a72f6c5 3003 for(Int_t iw = 0; iw < 14; iw++){
3004 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
649b825d 3005 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3006
3007 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 3008 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
649b825d 3009
3010 if(IsDataMC()){
3011
3012 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader(),0);
3013
3014 if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
3015 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) &&
3016 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3017 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3018 fhLambda0ForW0MC[iw][0]->Fill(energy,clus->GetM02());
1a72f6c5 3019 //fhLambda1ForW0MC[iw][0]->Fill(energy,clus->GetM20());
649b825d 3020 } // Pure Photon
3021 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3022 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3023 fhLambda0ForW0MC[iw][1]->Fill(energy,clus->GetM02());
1a72f6c5 3024 //fhLambda1ForW0MC[iw][1]->Fill(energy,clus->GetM20());
649b825d 3025 } // Electron
3026 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3027 fhLambda0ForW0MC[iw][2]->Fill(energy,clus->GetM02());
1a72f6c5 3028 //fhLambda1ForW0MC[iw][2]->Fill(energy,clus->GetM20());
649b825d 3029 } // Conversion
3030 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3031 fhLambda0ForW0MC[iw][3]->Fill(energy,clus->GetM02());
1a72f6c5 3032 //fhLambda1ForW0MC[iw][3]->Fill(energy,clus->GetM20());
649b825d 3033 }// Pi0
3034 else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3035 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3036 fhLambda0ForW0MC[iw][4]->Fill(energy,clus->GetM02());
1a72f6c5 3037 //fhLambda1ForW0MC[iw][4]->Fill(energy,clus->GetM20());
649b825d 3038 }// Hadron
3039
3040 }// Is MC
3041 } // w0 loop
3042
3043 // Set the original values back
3044 clus->SetM02(l0org);
3045 clus->SetM20(l1org);
3046 clus->SetDispersion(dorg);
3047
3048 }// EMCAL
3049
3050}
3051
3052
3053