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