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