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