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