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