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